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A Numerical Simulation of the Full 


Two-Dimensional Elect rothermal 
De-Icer Pad 

Konstanty C. Masiulaniec 


SUMMARY 


The ability to predict the time-temperature history of 
electrothermal de-icer pads is important in the subsequent design of 
improved and more efficient versions. These de-icer pads are installed 
near the surface of aircraft components, for the specific purpose of 
removing any accreted ice. The proposed numerical model can 
incorporate the full two-dimensional geometry through a section of a 
region (i.e., section of an airfoil, e:c.), that current one-dimensional 
numerical codes are unable to do. Thus, the effects of irregular 
layers, curvature, etc., can now be accounted for in the thermal 
transients. Each layer in the actual geometry is mapped via a 
body-fitted coordinate transformation nto uniform, rectangular 
computational grids. The relevant hea~ transfer equations are 
transformed and discretized. To model the phase change that might occur 
in any accreted ice, in an enthalpy formulation the phase change 
equations are likewise transformed and discretized. The code developed 
was tested against numerous classical and numerical solutions, as well 
as against experimental de-icing data on a UH1H rotor blade obtained 
from the NASA Lewis Research Center in Cleveland, Ohio. The excellent 
comparisons obtained show that this cole can be a useful tool in 
predicting the performance of current de-icer models, as well as in the 
designing of future models. 
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CHAPTER 1 INTRODUCTION 


The formation of ice on the exterior surfaces of aircraft has a 
considerable effect on flight performance, as it increases drag and 
decreases lift. Thus an aircraft must be designed with the equipment 
necessary for ice removal or prevention. Basically, aircraft ice 
protection systems can be classified as either anti-icing or de-icing. 

The anti-icing principle involves the prevention of ice formation 
on the protected area at all times. Typical anti-icing methods make use 
of chemicals and/or the passage of hot bleed air through channels below 
the surface on which ice formation is to be prevented. In contrast, 
de-icing involves the periodic removal of accreted ice by mechanical or 
thermal means. For ice removal systems, attention must also be given to 
uniform shedding of the ice. Itagaki [1 ] elaborates on the dangers of 
non-uniform shedding. Various de-icing methods that have been 
investigated include pneumatic boots and thermal techniques. The latter 
consists of cyclic heating of discrete elements by electrothermal means. 
The energy requirements are significantly less for de-icing systems than 
they are for anti-icing systems. From experimental studies, Stallabrass 
[2] concluded that the electrothermal method has the most advantages as 
a de-icing mechanism, although it does have some maintainability/ 
reliability problems. Werner [3] has also reported that the 
electrothermal de-icing technique is the most commonly used method, and 
that it has been applied to both fixed and rotary wing aircraft. 
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The objective of an electrothermal de-icing system is to raise the 
composite blade surface/ice interface temperature above the melting 
temperature of ice, resulting in a very thin interfacial layer of liquid 
which reduces the ice adhesion to the blade surface. Aerodynamic and/or 
centrifugal forces can then readily sweep the unmelted ice from the 
surface. A typical e lectrothermal de-icer pad is essentially a 
composite body consisting of (1) a metal substrate (the aircraft blade), 
(2) an inner layer of insulation, (3) a heating element, (4) an outer 
layer of insulation, and (5) an abrasion shield. Figure 1-1 depicts a 
two-dimensional cut-away view of the typical construction of an 
electrothermal de-icer pad, as well as a representative set of materials 
and thicknesses used for fabrication. The cross-section shown 
represents a view of the heater pad normal to the run of the heating 
e lements . 

The heating element usually employed in an electrothermal de-icer 
pad consists either of a woven mat of wires and glass fibers or of 
multiple strips of resistance ribbon. The gaps which exist between the 
heating elements can reduce the effectiveness of the heating pad 
de-icing performance, causing non-uniform melting of the ice. The two 
insulation layers, which usually consist of a resin impregnated ' glass 
cloth, serve to provide electrical insulation for the heating element. 

In order to direct more heat flow toward the ice layer, it is necessary 
to use a greater thickness for the inner insulation than for the outer 
insulation. The abrasion shield serves to protect the de-icer pad from 
rain erosion as well as dust/sand erosion, and to provide more uniform 
heating, thus minimizing cold spots above the heater gaps. 
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Layer 

Material 

Tiickness 

(Hr) 

Diffusivi ty 
(Ft 2 /I n) 

Substrate 

755-T6 Aluminum 

0.087 

1.65 

Inner 

Insulation 

Epoxy/Gl ass 

0.050 

0.0087 

Heater 

Ni chrome 

0.004 

0.138 

Outer 

Insulation 

Epoxy/Gl ass 

0.010 

0.0087 

Abrasion 

Shield 

304 Stainless 

0.012 

0.15 

Ice 


0.250 

0.0445 


Figure 1-1 Typical Materials and Construction of an Electrothermal 
De-Icer Pad 
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The ability to predict the performance of an electrothermal de-icer 
pad is essential to the design and subsequent fabrication of these 
units. To accomplish this, some method of determining the 
time-temperature history throughout the pad needs to be developed. 

Figure 1-2 provides a pictorial representation of an electrothermal 
heater section that is part of an airfoil, with some indication of the 
nature of the thermophysics involved. Clearly, the conduction of energy 
is three-dimensional in nature, and occurs in a curved, composite body. 
The temperature plot to the right of the figure provides a qualitative 
representation of a typical temperature distribution. The temperature 
is highest at the heater center, drops rapidly under the heater (where 
the insulation is thickest) and less rapidly in the direction of the ice 
(where the insulation is thinest). Development of an analytical model 
for such a problem is virtually impossible. A numerical model is more 
realizable, but even this is somewhat impractical, unless some 
simplifications are made to the geometry and the thermophysics. 

Figure 1-3 illustrates three alterations of the full de-icing problem, 
each having different degrees of problem simplification. The 
one-dimensional model is the simplest. In this model, all layers are 
assumed to be planes infinite in extent. The temperature at a given 
location is assumed to be constant throughout the plane containing that 
point. It is generally assumed that the layers are in perfect thermal 
contact and that they have constant material properties. 

Stallabrass [2] appears to have been the first to attempt a 
numerical solution of an electrothermal de-icing problem using a 
one-dimensional model. His numerical scheme used an explicit finite 
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difference method. Results agreed well with approximate analytical 
solutions for relatively short real times into the problem. To account 
for the effect of the phase change on the temperature transients within 
the composite blade, the node at the i ce-c^brasion shield interface was 
held at the freezing temperature until the estimated heat flux into the 
control volume containing the node was deemed sufficient to cause 
melting. 

Baliga [4] improved the numerical modelling of the same problem by 
handling the phase change heat transfer via a better approach, making 
use of the high heat capacity formulation. Marano [5] further improved 
upon Baliga 1 s numerical formulation by applying the so-called enthalpy 
method to model the phase change problem. Gent and Cansdale [6], solved 
the same problem for conduction only (no phase change), and obtained 
nearly the same results as Marano for concuction only. 

The two-dimensional problem, represented by the middle schematic in 
Fig. 1-3, was solved by Chao [7] and DeWitt, et al. [8]. Chao's work 
was a direct extension of Marano's one-dimensional numerical formulation 
and procedures to two dimensions. Of fundamental importance, the effect 
of the heater gap width on de-icing performance was studied 
numerically • 

Leffel [9] provided detailed experimental results of the thermal 
transients induced by an electrothermal de-icing unit on a UH1H 
helicopter rotor blade section. These experimental results were used to 
validate the codes developed by Chao and Marano. The experimental 
results revealed that when the layers of e helicopter blade are 
sufficiently thin, and the curvature sufficiently gradual, Marano's 
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one-dimensional code yields excellent results over most of the blade. 
Furthermore, it was found that there are two regions of potentially 
substantial inaccuracies (depending on heater wattages, material 
properties, etc.). These are at the immediate edges of the heater 
banks, and in the region of large curvature at the leading edge of the 
blade that wraps around the nose block. 

Chao*s code can model the heater edges, but it can not handle the 
variable thickness introduced by the nose block, nor the high degree of 
curvature. Thus, it is necessary to develop a model that can account 
for these difficulties. This development is pictorially represented by 
the third schematic of Fig. 1-3. The creation of a numerical code 
capable of accurately predicting electrothermal de-icer pad thermal 
transients in more complex regions of the blade is the topic of this 
work • 

Overview of Strategy to Solve Problem 

The primary modelling difficulty that must be faced in this study 
is that due to the irregularity of the blade-layer geometry. There are 
essentially three approaches that can be taken to account for irregular 
body curvature. The first is to overlay the irregular geometry with a 
regular grid. Those points not falling directly on the boundary will 
require an interpolating scheme that must be incorporated either 
directly or indirectly into the computational algorithms of the field 
equations. There will undoubtedly be some degree of inaccuracy 
introduced into the solution as a consequence. If the solution desired 
is in a region reasonably far removed from the boundary, this approach 
will generally yield excellent results. If, however, the solution 
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desired is located at the boundary, the inaccuracies introduced may be 
unacceptable. Such is the case in this problem, where the critical 
temperatures and the initial change of phase occur at the boundary 
formed by the abrasion shield and the ice layer. Moreover, this 
approach is computationally quite time consuming. 

The second approach is the finite element method. This technique 
has become fairly common and affords many advantages in problems having 
irregular geometries. The standard method needed to achieve a solution, 
however, is by an inversion of the appropriate system matrix. For 
problems requiring a large number of elements, a large matrix develops. 
The inversion of this matrix is generally performed by iterative means, 
thus affording no computational advantage over other techniques. Also, 
the formuation of the problem required for more conplex governing 
equations will involve the inversion of a multiple number of large 
matrices for a "matrix statement" equivalent [10], [11] of the field. 

In recent years, a finite difference alternative has arisen that 
accurately models a boundary of irregular shape. By this method the 
body (or bodies) in the physical domain is numerically transformed into 
a rectangular region in the computational plane. For a layered body, 
the transformation produces a rectangle having the same number of layers 
as the composite in the physical plane. This mapping procedure, known 
as body fitted coordinate generation, was first developed by Thompson, 
et al [12], [13], [14], [15], and has become widely used in numerical 

simulations of field problems [16], [17], [18], [19], [20]. The 

principal advantage of the procedure is that any set of equations may be 
numerically solved in a rectangular region on a uniformally spaced grid 



system. Thus, numerical interpolation between any irregular boundary 
and adjacent interior grid points can be avoided. Figure 1-4 depicts 
the mapping strategy involved for the problem at hand. The numerical/ 
computational strategy needed essentially reduces to a finite difference 
solution of a series of stacked, rectangular slabs. The primary 
disadvantage of this technique is the necessity of spacially 
transforming those portions of any relevant equations having a spacial 
dependency. Depending on the equations to be transformed, the resulting 
set of equations may become much more complex. Obviously, this has the 
potential of making a numerically stable set of algorithms more 
difficult to obtain. 

This work represents the first known attempt at solving a layered 
heat transfer problem that includes phase change in an irregular 
geometry. The following chapter develops the spacial transformation 
equations, and the operators needed to transform the governing 
equations. The conduction equations are then transformed in Chapter 3, 
with the phase change equations being transformed in Chapter 4. These 
transformed equations are then used to simulate an iced airfoil, and the 
predictions are compared with experimental test data. 
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Figure 1-4 Pictorial Representation of how the 
Composite Blade Geometry is Mapped 
with Body Fitted Coordinate Transform 
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CHAPTER 2 - BODY-FITTED TRANSFORMATION 

In order to obtain a grid in the transformed plane a method with a 
system of generating equations needs to be developed. The system of 
equations that was used to generate the grid was Laplace's equations. 
Virtually any partial differential equation can be used. Thompson, et 
al. [15] presents alternatives by adding additional terms to Laplace's 
equation that skews the grid in a desired direction where exceptionally 
large gradients of a variable are expected. Laplace's equation was 
chosen since this equation type is closest to the form of the equations 
that was solved in the transformed plane. Thus, there are two 
rectangular coordinate systems that are interrelated through the 
equations that generate the grid in the transformed plane: x-y in the 

real plane, and E-n in the transformed plane. 

It should be noted that for composite bodies, each region must be 
mapped separately to insure accurate modeling of the true geometry. The 
initial step of the mapping is to assign boundary points in the real 
plane, which become boundary conditions in the transformed plane. Thus 
by mapping each region separately, more accurate geometric modeling is 
achieved. The assigned boundary points for two regions having a 
boundary in common must be the same to insure a continuous grid 
between regions. This makes the subsequent application of boundary 
condition equations much less complex. 
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Since the generating equation is initially written in the real 
plane, an exchange of variables must occur to obtain an equivalent 
set of equations in the tranformed plane governing the distribution of x 
and y coordinates in that region. This is the set of equations that is 
subsequently used to generate a grid. These equations will be developed 
first, followed by a physical interpretation of the transformation. 


Derivation of Body-Fi tted Coordinate Transform Equations 

The partial differential equation used to generate the grid in the 
transformed plane, Laplace's equation, is written as 

2 2 
S£ 9£ 

— 3 +—3=0 (2-1 ) 

ay 


.2 2 

a r\ an 

--- + — - = 0 

a x z ay^ 


( 2 - 2 ) 


The first step in performing the required exchange of variables is to 
establish derivatives for the inverse trarsform. In general notation 
the inverse transform is represented as 

x = f <C f n) (2-3) 

( 2 - 4 ) 

(2-5) 


y = g ( C, n) 

The total derivative of each variable is 

9f 3f 

dx = -- d £ + — d ti 

H 



( 2 - 6 ) 
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3g 3g 

dy » -- d 5 + — d n 

3 5 3n 

Multiplying Eq. (2-5) by 3g/3n and Eq. (2-6) by 3f/3r) and 
subtracting yields 


3g 3f / 3f 3g 3g 3f \ 

— dx dy = ] d5 (2-7) 

3n 3n \3C 3n 35 9n/ 


Similarly, Eqs . (2-5) and (2-6) may be multiplied by 3g/35 and 3f/9£, 
respectively, and subtracted to give 

3g 3f / 3g 3f 3f 9g\ 

-- dx - — dy = ( — -- - -- — ] dn (2-8) 

35 35 3n 3? 3n/ 


The Jacobian of the transformation may be defined as 


3f 3f 

3 5 3 n 3f 3g 3g 3f 

J = (2-9) 

3g 3g 3 5 3n 35 3n 

3 5 


Inserting this into Eqs. (2-7) and (2-8) and rearranging produces 


1 3 g 1 3f 

d 5 = dx dy (2-10) 

j 3n J 9n 


1 3g 1 3f 

d n = - - -- dx + - -- dy (2-11) 

J 35 J 35 


Since 5 and n are transforms of x and y, the derivatives of 5 and n are 


written as 



differential operators, i.e., 
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If the following are defined, 
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Equation (2-22) becomes 


9 ^ 3 ^ 1 / 9 1 ^ 9 ^ \ 

___ + ___ = __ (a — - 26 + Y — - ) 
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(2-26) 


Equation (2-26) can now be used directly to obtain the spacial field 
equations needed to generate the coordinate lines in the physical 
plane. Beginning with Laplace's equation in the real plane as the 

generating equation, these equations are 

2 2 2 

9 x 9 x 9 x 

a — - -26 + Y — ^ = 0 (2-27) 

3£ 3£3n 9n 

2 2 2 
3 y 3 y 3 y 

a -2 6 + Y — 2 - 0 (2-28) 

H 3£9n 9n 

Because Eqs . (2-27) and (2-28) have the variable coefficients alpha, 

beta and gamma in common (which are derivative functions of x and y), 


they are coupled and must therefore be solved simultaneously. 

In order to formulate an equality of the normal flux boundary 
condition in the transformed plane, a vector derivative transformation 
needs to be developed. The unit normal to any graph f(x,y) = c is given 
by 

♦ w 


n f = 7" 7 

I I 

Consequently the unit normal to n and £ coordinate lines are 


(2-29) 



The solution of Eqs. (2-27) and (2-28) results in a distribution of x 
and a distribution of y in the solution plane. Thus, in order to 
represent a unit normal vector in this plane, its 'orientation* in the 
real plane must be accounted for. The general operator in the real 
plane for a directional derivative is 

3 3 4 . 

V = — i + — j (2-32) 

3x 3y 

Applying Eqs. (2-18) and (2-19) to Eq . (2-32), the general operator 
in the transformed plane for a directional derivative is 




(2-33) 


It should be noted that the directional derivatives are needed in the 
real plane, but are being represented by information that is contained 
in the transformed plane. Thus, the symbols i and j refer to the real 
plane and do not need to be transformed. 


Using Eq. (2-33), Vn and clearly become 

1 / 3y + 3x -A 

Vn = - [ i + — j J (2-34) 
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In the next section a physical interpretation of the coefficients alpha, 
Eq. (2-23), and gamma, Eq. (2-25), and how they are related to arc 


lengths of cells in the physical plane will be presented. It will be 



shown that 
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Hence, utilizing Eqs. (2-34) through (2-37), Eqs. (2-30) and (2-31) 
become 
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The temperature gradient at any point in the point plane is defined by 

Vt = — i + -- j (2-40) 

3x 9y 


Applying Eqs. (2-18) and (2-19) yields the temperature gradient along 
any point in the transformed plane as 
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(2-41) 
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The gradient of temperature normal to a line of constant n becomes 

3t -► 



1 

j 



(2 


The gradient of temperature normal to a line of constant £ becomes 

9t + 

= n { • Vt 



- 42 ) 


( 2 - 43 ) 
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Equations (2-42) and (2-43) are those needed in formulating the normal 
heat flux condition f or a boundary common between two regions. 

Physical Interpretation and Use of the Body Fitted Coordinate 
Transf orma tion 

Although field equations are written in terms of physical plane 
coordinates (x,y; r, A; etc.), their analytical or numerical solution is 
generally only convenient for a limited number of smooth, regular 
geometries. If the geometry is irregular the solution in the physical 
plane can become quite tedious. An irregular body can be mapped into a 
much simpler shape by the use of a body fitred coordinate transfor- 
mation, as shown in Figs. 2-1 through 2-3. 

If the composite body is simply connected it will have an interior 
solid which requires a slightly different treatment. This is also 
illustrated in the bottom two rectangles of Fig. 2-3. Care must be 
taken to match the edges of the rectangle w.th corresponding portions of 
the interface for the layer which surrounds the solid. 

The overall result of the transformation is the creation of a 
series of connected rectangular strips. To generate the coordinate 
system within these strips, two elliptic partial differential equations, 
Eqs. (2-1) and (2-2), with Dirichlet boundary conditions are solved. 

The field equations used to develop the coordinate system are Eqs. (2-27) 
and (2-28), which are the end result of the exchange of variables needed 
for the transformed plane. The variable coof f icients alpha, beta, and 
gamma are defined by Eqs. (2-23) through (2*25). 

Figures 2-4 through 2-6 conceptually s iow the procedure needed to 
generate a grid, and a sample of what the final results might yield. 



22 



Figure 2-1: 


A General, Irregularly Shaped Composite 
Body . 




Figure 2-2: Division and Separation of an Irregularly 

Shaped Composite Body for Subsequent 
Mapping in the Transformed Plane. 
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A set of x-y coordinate pairs need to be identified for the boundaries 
of a given region in the physical plane. These values are then assigned 
to boundary points for a rectangular region in the transformed plane. 

The distribution of x,y coordinates in the £, h plane are tied together 
by the partial differential equation used to generate the grid. 

Equations (2-27) and (2-28) are then solved simultaneously on the field, 
with x and y coordinate values as the unknowns. Upon convergence, the 
solution yields a distribution of x and y coordinates within the field, 
thereby accomplishing the goal of generating a grid. If the resulting 
solution is in turn mapped onto the physical plane, its appearance would 
take on the shape of a collection of quadrilaterals that approximate the 
original regions. Thus, it can be seen that an improved approximation 
is achieved as the number of boundary points chosen is measured, 
especially on those portions of the boundary that have a high degree of 
curvature. The numerical solution requires the same number of arbitrary 
points on the inner boundary as on the outer boundary to maintain a 
rectangular grid in the transformed plane. Also, composite or layered 
bodies should have the same boundary points for shared boundaries, as 
otherwise a temperature solution involving more than one region would be 
very difficult to numerically implement. Figure 2-6 shows that for a 
conduction heat transfer problem in the transformed plane there are 
always three solution fields. The x and y solution fields are 
determined first (which generates the grid), and then on this grid a 
solution is obtained for a field equation(s) of interest. Thus, each 
nodal coordinate for a specified value of £ and n has at least three 
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oundary Values (Dirichlet Boundary Condition) in the Transformed Plane. 






Figure 2-6: The Three Solution Fields Present for each Region. 
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values associated with it: an x coordinate, a y coordinate, and, in the 

case shown, a temperature* All of these values at each node must then 
be used to obtain temperature mappings in the physical plane* 

The coefficients alpha, beta and gamma are functions of x and y* 

In the final solution to the grid generation, there is also a 
distribution of these three coefficients throughout the field. Thus, 
these coefficients can change in value from point to point. The 
physical significance of all the transf orma tion coefficients (including 
the Jacobian) is illustrated in Figs. 2-7 through 2-10. Figure 2-7 
shows an arbitrary quadrilateral drawn onto a square grid in the 
physical plane. The midpoints of the opposite sides of the edges are 
connected, which give essentially an average "length" and "width" of the 
quadrilateral. These dimensions are lableo "b" and "a" in Figs. 2-8 and 
2-9, respectively. All of the derivatives of x and y with respect to £ 
and n are taken as discrete, constant values through a quadrilateral 
along the lengths a and b. It is inherently assumed that the lengths a 
and b correspond to lines of constant K and n from the transformed 
plane. The intersection of lines a and b represent, then, the 
intersection of the lines £ and n for a given point that was part of the 
iterated field in the transformed plane. 

From reference to Figs. 2-7 through 2-10, it can be seen that the 
Jacobian of the transformation is the area of the quadrilateral in the 
physical plane. The square root of alpha Isa measure of the average 
quadrilateral length along the n direction in the physical plane. The 
square root of gamma is a measure of the average quadrilateral length 
along the £ direction in the physical plane. And, finally, beta is a 




, Figure 2-7: Spacial Relationship of a Computational Cell in 

the Transformed Plane Mapped onto the Real Plane. 




Figure 2-8: Physical Significance of the Transformation 

Parameter 'Alpha' in the Real Plane. 
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Figure 2-9: Physical Significance of the Transformation 

Parameter 'Gamma' in the Real Plane. 


33 



= -.447213 


= ab coss 


/ 40( 200) 


Figure 2-10: Physical Significance of the Transformation 

Parameter 'Beta' in the Real Plane. 
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measure of the non-orthogonality of the quadrilateral as measured by the 
mapping of the £ and n lines on the corresponding cell/quadrilateral in 
the physical plane. Thus, the essential dimensions describing the 
approximate shape of the quadrilaterals in the physical plane (which are 
the square, uniform computational cells in the transformed plane) are 
contained within the transformation parameters alpha, beta, gamma and 
Jacobian. 

Numerical Modelling of Body-Fit Transformation Equations 

Solving the transformation equations in the computational plane is 
essentially the same as solving any two field equations that are coupled 
with Dirichlet boundary conditions. Equations (2-44)* through (2-50)* 
show how the two field equations have been discretized. All second 
order derivatives are approximated using three point central 
differences, and all first order derivatives using two point central 
differences. Note that the coupling of the equations occurs through the 
transformation coefficients alpha, beta, and gamma. Equations (2-45) 
and (2-47) can be solved for x(i,j) and y(i,j) to obtain the required 
algorithms. Figure 2-11 illustrates a representative portion of the 
computational grid that Eqs . (2-44) through (2-50) were solved on, 
explaining the indices i and j. 

The iteration strategy is to first sweep through the x field, 
holding all values of y at the value of the most recent iteration, and 
then sweep through the y field, holding all values of x fixed from the 
most recent iteration. The values of the transformation coefficients 


*Because of their length, these equations are shown on the following 
page. 
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r\ ^ ^ ^ 

OX ox c x 

ct -2 8 + Y --- = 0 (2-44) 

9? 2 9£9n £t/ 


X i+1,j+1 - X i-1,j+1 X i+1,j-1 " 


a(x. , . - 2x. . + x_. , , ) -2 fi 


i+1 * j 


. + A . . 

,9 1-1*9 


+ Y (x. . , - 2x . . + 

1*9 + 1 1*9 


9 2 y 


H 


a(v. . - 2y. . + y. , 

J i + 1, 9 1*9 1-1* 


) -28 


2 2 

2 


-i> = 0 

(2-45) 

9 2 y 9 2 y 

. 28 + '• — - = 0 

9E;9n 3n 

(2-46) 

y i+1 , j+1 " y i-1,j+1 y i+1 , j -1 " 

y i - 1 * j - 1 

2: 2 


+ Y (y. . , - 2y. . 

1,9 + 1 1*9 


y. . J = 0 

1 * 9-1 


(2-47) 


a = 




z 

(2-48) 


3 y 3 y 3 x 3 x 

8 n 3 £ 9 n 3 ^ 


y i,j+1 ' Y i , j-1 y i+1,j ~ Y i-1*j 


2 


2 


x i + 1 * j “ x i~ 1 * j 

+ 

2 


y i+1,j " y i-1*j 


2 


( 2 - 49 ) 



(2-50) 
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Figure 2-11: Grid Index Notation for Solution of Spacial 

Transform Algorithms in Transformed Plane 
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are updated for each point for every iteration to reflect the most 
recent values of x and y, as the sweep through the field progresses. 

This procedure continues until all values cf x and y are settled i.e., 
they no longer change to within some specified error. 

Prior to the beginning of the iteration, a linear distribution 
between any two opposing boundary values wc ; s assigned to all the field 
points as an initial guess. This not only reduces the needed CPU time, 
but also assures that there will be no divide check errors during the 
first sweep of the field. Also, once convergence has been achieved, the 
values of the transformation parameters for the boundary points are 
calculated for use in the boundary condition algorithms. The derivative 
values needed are calculated using one-sided differencing. Appendix A 
contains a computer listing of the key subroutine written in Fortran. 

Verification of the Transformation Code/Algorithms 

In order to verify the accuracy of the coordinate mapping, use will 
be made of the fact that an analogy exists between the spacial transform 
solution and a solution of a field equation that has the same form as 
the mapping equations. To make the comparison, Laplace* s equation for 
steady state conduction was used: 

3 Z T 3 

+ — = 0 ( 2 - 51 ) 

3x 1 3y ^ 

This equation was solved analytically for the case of two concentric 
cylinders having fixed temperatures at the inner and outer surfaces. 
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The solution to the problem is well known, and is given in terms of 
radial coordinates by 


r 

In — 
r o 

T(r) = (Ti-T 0 ) + T 0 (2-52) 

r i 
in — 

r o 

Using this equation the radial location for a specified temperature may 
be written 


T-T 

o 



For an inner and outer radius value of 1 and 9 feet, respectively, and 
an inner and outer temperature of 400 and 0 degrees Farenheit, 
respectively, the locations of the isotherms in 40 degree Farenheit 
increments were calculated, and are listed in Table 2-1. Those radius 
values for temperature increments of 10 degrees Farenheit are listed in 
Table 2-2. From a numerical point of view, these two tables represent 
exact values for a fine and coarse mesh. Listed along side these values 
are the final results of a spacial transform to map the given circular 
cylinder. The output for the mapping on a 41 x 41 grid can be found in 
Appendix B, and the output for the 11x11 grid in Appendix C. Note 
that each coordinate point has an x and y value associated with it. 
Knowing that the boundary points were assigned with the region centered 
on the origin of an x-y coordinate system, these values can be easily 
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converted into an equivalent radial location* This task is much simpler 
if values are chosen along the x or y axis. The radial locations can 
then be picked off directly. 

Because the generating equations used are Laplacian in form, there 
is an exact analogy between the resulting and y distributions and the 
location of isotherms, where those isother ns are equally spaced (i.e., 
each grid point outward from the inner boundary represents an equal drop 
in the number of degrees Farenheit). The percent error between the two 
solutions presented in Tables 2-1 and 2-2 clearly shows increased 
accuracy with a finer mesh. The error for an 11 x 11 grid drops from a 
maximum of 7.13% to a maximum of 0.40% for a 41 x 41 grid. The error, 
as would be anticipated, increased with increased distance from the 
boundary where the temperature is known. 
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TABLE 2-1 

ACCURACY OF ITERATED SPACIAL VALUES 

FOR TWO 


CONCENTRIC 

CIRCLES USING A 



41 x 

41 GRID 


Isotherm 




Value 

Ana lysis 

Program 

% 

( °F ) 

Radius Value 

Radius Value 

Error 

0 

9.0000 

9.0000 

0.0000 

10 

8.5189 

8.51 55 

0.0399 

20 

8.0636 

8.0573 

0.0781 

30 

7.6326 

7.6240 

0.1126 

40 

7.2246 

7.2141 

0.1 453 

50 

6.8385 

6.8263 

0.1 784 

60 

6.4730 

6.4596 

0.2070 

70 

6.1 270 

6.1 1 27 

0.2333 

80 

5.7995 

5.7845 

0.2586 

90 

5.4895 

5.4741 

0.2805 

1 00 

5.1 961 

5.1 804 

0.3021 

1 1 0 

4.91 84 

4.9026 

0.3212 

1 20 

4.6555 

4.6397 

0.3393 

130 

4.4067 

4.391 1 

0.3540 

1 40 

4.1711 

4.1 558 

0.3668 

1 50 

3.9482 

3.9333 

0.3773 

1 60 

3.7371 

3.7227 

0.3853 

1 70 

3.5374 

3.5235 

0.3929 

180 

3.3483 

3.3350 

0.3972 

1 90 

3.1 694 

3.1 566 

0.4038 

200 

3.0000 

2.9878 

0.4066 

21 0 

2.8396 

2.8282 

0.401 4 

220 

2.6878 

2.6771 

0.3980 

230 

2.5442 

2.5342 

0.3930 

240 

2.4082 

2.3989 

0.3861 

250 

2.2795 

2.2709 

0.3772 

260 

2.1 576 

2.1 497 

0.3661 

270 

2.0423 

2.0351 

0.3525 

280 

1 .9331 

1 .9266 

0.3362 

290 

1 .8298 

1 .8240 

0.3169 

300 

1 .7320 

1 .7268 

0.3002 

310 

1 .6394 

1 .6348 

0.2805 

320 

1 .551 8 

1 .5478 

0.2577 

330 

1.4689 

1 .4655 

0.2314 

340 

1.3903 

1.3875 

0.201 3 

350 

1.3160 

1.3137 

0.1 747 

360 

1.2457 

1.2439 

0.1 444 

370 

1.1 791 

1.1 778 

0.1 102 

380 

1 .1 1 61 

1 .1 153 

0.0716 

390 

1.0564 

1 .0560 

0.0378 

400 

—ft 

• 

o 

o 

o 

o 

1 .0000 

0.0000 
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TABLE 2-2 

ACCURACY OF ITERATED SPACIAL VALUES 
CONCENTRIC CIRCLES USING AN 11 x 11 

FOR TWO 
GRID 

Isotherm 
Value 
( °F) 

Analysis 
Radius Value 

Program 
Radius Value 

% 

Error 

0 

9.0000 

9.0000 

0.0000 

40 

7.2246 

7.0175 

2.8665 

80 

5.7995 

5.51 37 

4.9280 

1 20 

4.6555 

4.3631 

6.2807 

1 60 

3*7371 

3.4757 

6.9947 

200 

3.0000 

2.7860 

7.1 333 

240 

2.4082 

2.2463 

6.7228 

280 

1 .9331 

1 .821 1 

5.7938 

320 

1 .551 8 

1 .4841 

4.3626 

360 

1.2457 

1.2154 

2.4323 

400 

1 .0000 

1 .0000 

0.0000 


42 


Since the spacial transform essentially maps an irregular region 
into a collection of quadrilaterals, the distribution of transformation 
coefficients can be used as a check on the accuracy of modelling the 
entire geometry* Table 2-3 shows this comparison, where the overall 
values of thickness, circumference and area of the annulus are 
calculated by summing the appropriate coefficient values* The largest 
error for both grid sizes is in the approximation of the circumference. 
This is to be expected, since the curvature is being approximated by a 
series of straight lines. Again, the finer mesh provides a much more 
accurate modelling of the geometry. The largest error is reduced from 
6.46% to 0,419%. The error in radius for a 41 x 41 grid is slightly 
larger than for an 1 1 x 1 1 grid. This is attributable to round-off 
error, as only single precision was used in the numerical code. 

This comparison and verification demonstrates the accuracy of the 
spacial transformation subroutine. However, it also points out the 
possibility of large errors entering the solution of a field equation if 
the grid that is used is not sufficiently fine to accurately model the 


region 
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TABLE 2-3 ACCURACY OF SPACIAL TRANSFORM BASED ON PHYSICAL 

SIGNIFICANCE OF TRANSFORMATION COEFFICIENTS COMPARED 
TO ACTUAL VALUES OF THE REGION IN THE PHYSICAL PLANE 


Physical 

Significance 


Radius 

(R 0 - R i> 


Circumf erence 
(2 ttR 0 ) 


Area 

2 2 

( TT[R 0 _R i ] ) 


Parameter 


I /oT 


Actual 

Value 


T 


Actual 

Value 


I J 


Actual 
Va lue 


11x11 GRID 


7.9999 


8.0000 


52.9006 


56.5486 


235.1 146 


251.3274 


0 . 001 % 

Error 


6.451% 

Error 


6.450% 

Error 


41 x 41 Grid 


7.9880 


8.0000 


56.3162 


56.5486 


250.2978 


251.3274 


0.025% 

Error 


0.419% 

Error 


0.409% 

Error 













CHAPTER 3 TRANSFORMATION OF CONDUCTION PROBLEM 


Transformation of Governing Conduction Equations 

Any field equation solved in the transformed plane containing terms 
having a spacial dependency must have those terms transformed in the 
same fashion as the spacial coordinates £, n were developed. The form of 
the conduction energy equation used in this study is: 

3t 3 / 9t\ 3 / 3t\ 

0Cp k ( k — J + q ( 3-1 ) 

3t 3x\^ 3x J 3y ^ 3y J 

In the problems considered the thermal properties of the materials are 
assumed to be constant; hence, Eq . (3-1) can be rewritten as 

3t / 3 * T 3 ^T 

0Cp ' k + i ; 5 

To transform the above equation, the transformation operator developed 
in the previous chapter, Eq. (2-20), is required, i.e., 


+ q 


(3-2) 


3x 3y 


1 


z 

J 


a — - -2R + 




3£3n 



(3-3) 


where 
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3x 3x 3y ty 

ft + — — (3-5) 

9S 3n ar 3n 



3x 3y 3y ;bc 

and J = -- — - — (3-7) 

as 3n as 


Applying Eq . (3-3) to Eq . (3-2) yields: 


3t 


3*T 


3^T 


3 ' 


Pc r 


ex --- -2 ft + y —3 I + q 

as 3S3n 3n ' 


(3-8) 


3 t J 

This is the transformed energy equation that must be solved 
in the transformed plane. 


Transformation of Boundary Condition Equations 

There are two basic methods that can be used for treating an 
interfacial boundary condition. One is to equate normal heat fluxes and 
the other is to perform an energy balance on a control volume that 
includes the boundary, i.e., 



VT*n)ds + 



(3-9) 


For a square of rectangular geometry, application of Eq. (3-9) is 
relatively straightforward. An approximation of Eq. (3-9) can be 
derived without any special treatment of tie interface, and the 
resulting discretized equation can be cast into simple algebraic form. 



46 


To apply the above energy balance to a square cell or control 
volume in the transformed plane is not as straightforward, however. 
Equation (3-9) must undergo the same spacial transformation as 
described above. To apply the resulting transformed equation to a 
control volume in the computational plane is possible, but the added 
complexity makes the task quite prohibitive. 

Another way of "visualizing" the task is to apply Eq. (3-9) 
directly to a control volume in the physical plane that has been 
"unmapped" from the computational plane. Figure 3-1 shows how the 
resulting cell might appear in the real plane. As can be seen, each 
cell around the boundary has a different configuration. 

Computationally , customizing Eq. (3-9) to each cell around the boundary 
is not easily accomplished, even with a knowledge of the approximate 
cell geometry and configuration from the transformation parameters at 
each point. Accuracy needs to be preserved to as great an extent as 
possible, while maintaining a reasonable degree of computational 
simplicity. Thus, Eqs. (2-42) and (2-43), which can be used to equate 
normal heat fluxes through a given boundary point along a line of 
constant £ or n, shall be used for the boundary conditions between two 
adjacent regions. 

The boundary conditions for two regions in perfect thermal contact 
require that there be equality of temperature and normal heat flux: 


T 


1 


I 



(3-10) 


*1 




(3-11 ) 
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On the other hand for an exposed surface, energy conducted to the 
surface is convected away; therefore, 


9t 

-k — = h (T _ T J 

s 


(3-12) 


Since the boundary equations contain normal derivatives that are 
functions of space, a transformation operator is required that will 
provide an expression for the temperature derivative normal to the 
surface at a specified point in the real plane. From Chapter 2, these 
are 


9 

1 

/ 3t 

*r\ 


Sn 

To a line of 

( Y 7, 

”«) 

(3-13) 


constant h 


a 


3n 


To a line of 
constant K 


1 

J /aT 



3t 3t\ 



(3-14) 


On the presumption that the boundary is along the upper or lower edge of 
the transformed plane (a line of constant n), Eqs . (3-11) and (3-12) 
become , respectively , 






(3-15) 



P &r 

-j 

as 



h(T s - T J 


-k 


(3-16) 
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Equations (3-15) and (3-16) are the transformed layer surface boundary 
equations that must be solved in the transformed plane* 

Numerical Differencing of Conducting Equat ions 

Figure 3-2 shows a section of the grid in the transformed plane* 

In order to accurately difference the cross-derivative terms in the 
diffusion equation, all eight nodes surrounding the node being iterated 
for must be involved in the differencing* Applying second order central 
differences to all derivatives in Eq. (3-9) yields 



kAt 


k+1 

(T i,j 



(1 -M) 



(T i+i,j 


k k 

2T i , j + T i-1 , j ^ 


p. . . 

T i+ 1 ( j +1 

2 


k k 

T i+1 , j-1 + T i-1 , j-1 




+ Y i,j (T i,j + 1 “ 2T i / j + 


k 

r i,j-1 5 


+ M 



k+1 

(T i+1 , j 


2T. 


k+1 
if j 


+ T 


k+1 
i - 1 / D 


) 


R i,j k+1 k+1 <+1 

(T i+ i,j+i “ T i+1/j-1 + T l- 1 F j -1 

2 


k + 1 

T i-1, j+1 


J 


+ Y. 


k+1 


k+1 k+1 

2T i,j + T i, j-1 ’ 


+ 


2 

J i,j « 


k 


(3-17) 


It should be noted that a weighting factor M has been used in order that 
the algorithm may become anything from purely implicit (M equal to zero) 
to purely explicit (M equal to one). When M is equal to 0.5, the 
algorithm becomes the Crank-Nicholson scheme. Appendix D contains a 



50 














51 


listing of the subroutine that solves for the conduction problem 


utilizing Eq. (3-17). 

Upon differencing, the boundary condition for equality of normal 
heat flux at a solid boundary becomes: 

T i+1 , 1 ~ T i-1 , 1 \ 

r / 

upper 


y. 

1 , s 


v i,s 




< T i ,2 - * 1,1 ) ~ 




^i , j max 



y. 

l , j max 


JY. ! 

,jmax i,jmax 


( T i 


, j ma ;< 


T i , j max- 1 


) 


, j max 


J i , j max i , j max 


T i+1,jmax " ^i-1,jmax 
2 


(3-18) 


lower 


The algorithm for a convective boundary condition is similar to 

Eq. (3-i8) # except that the left or right hand side is replaced by 

h (T. - T m ) or h CT ; " T J, depending on whether the convective 

1,1 l / j max 

boundary is on the upper or lower surface c-f the grid. Appendix E 
contains a listing of the subroutine that solves for the normal flux 
boundary condition using Eq . (3-18). 

Verification of Numerical Modelling 
A. One Zone, Steady State 

The problem used to verify the algorithm for conduction in a single 
layer under steady state conditions is conduction in offset cylinders. 
The exact solution to this problem is found, in Eckert and Drake [21]. 
Computed isotherms for a cylinder with an outer radius of 0.3048 meters 
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(1.0 feet) held at 37.78 degrees C (100 degrees F), an inner radius of 
0.0457 meters (0.15 feet) held at 537.78 degrees C (1 000 degrees F), and 
an eccentricity of 22.86 millimeters (0.75) feet are shown in Fig. 3-3 
(straight line). For a grid consisting of 17 x 17 nodes, the numerical 
values are in excellent agreement with the analytical values. 

The largest error to be found between any two values was less than 0.1%. 
Thirty— seven seconds of computing time was needed to obtain a converged 
result on an IBM 4341. 

The eccentric cylinder problem was also solved analytically for the 
case involving uniform internal heat generation by El-Saden [22]. This 
problem was numerically solved for the case where energy was being 
uniformly generated at the rate of 20.7 kilojoules per second - cubic 
meter (2000 BTU per hour - cubic foot). Figure 3-3 compares the exact 
to predicted temperatures on the lines £ = 1, £ = 6, and £ = 8. These 
coordinate lines represent cuts in the transformed plane at 0, v/2, and 

tt respectively. Again there is excellent agreement, with the maximum 
error being less than 0.5 %. 

B. Multiple Zones, Steady State 

The steady state temperature distribution in four concentric 
cylinders was obtained and compared to the elementary analytical 
solution for this problem. Table 3-1 lists the parameters used for this 
problem. 

Table 3-2 presents the predicted temperatures as compared to the 
analytical values at three interfaces for various convergence criteria. 
The surface temperature of the inner cylinder was held at 1000 degrees 
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With and Without Heat Generation 



54 


i 

Inner 

Radius 

Outer 

Radius 

Thermal 

Conductivity 

Layer 

1 

Rj_(TTim) 

R Q ( mm ) 

k(kJ/hr-m-°C ) 

1 

1 000 

800 

155.77 

2 

800 

700 

249.23 

3 

700 

500 

18.69 

4 ! 

500 

300 

93.46 


Table 3-1 : Parameters Used in the Multiple Zone Steady 

State Verification Problem 




Table 3-2: Code Verification Results - Interface Temperature 

Solutions for Concentric Cylinder Problem 


Analytical 

Num. 1 

Num. 2 

Num. 3 

Num. 4 

Num. 5 

( °c) 

(°c) 

( °c) 

(°C) 

(°c) 

( °c) 

1 00,00 

100.00 

k 

O 

O 

• 

o 

o 

100.00 

1 00.00 

100.00 

1 50,69 

1 73.25 

158.3 

1 55.6 

154.0 

1 53.3 

1 69.65 

196.0 

1 78.5 

1 75.4 

173.5 

1 72.74 

806.60 

81 5. 1 

809. 72 

808.9 

808.4 

808.1 5 

1 000.00 

1 000.00 

1 000.00 

1 000.00 

1000.00 

1 000.00 

m 

ii 

it 

ii 

ii 

ii 

ii 

ii 

it 

ii 

n 

ii 

it 

m 

ii 

ii 

n 

n 

it 

ii 

n 

ii 

n 

ii 






Spatial 

Convergence 

0.0001 

0.0001 

0.0001 

0.0001 

0.0001 

Ttemp. Field 
Convergence 

0.1 

0.1 

0.1 

0.1 

0.1 

Boundary 

Convergence 

0.1 

0.1 

0.05 

0.02 

0.01 

CPU Time 
IBM 4341 

1 : 29 

3:31 

3:52 

3:58 

4:00 
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C, while the surface temperature of the outer cylinder was held at 100 
degrees C. 

The error for the innermost boundary is 2.61 degrees C, for the 
middle boundary 3.09 degrees C, and for the outermost boundary 1.55 
degrees C. The largest error for numerical solution 5 is less than 1.8 
percent. The above solutions were obtained on an 11 x 11 grid for each 
layer. Figure 3-4 displays a graphical representation of the numerical 
solution. In light of the relatively coarse grid used, the agreement is 
excellent. 

C. One Zone, Time Dependent 

Jakob [23] developed exact solutions for three one-dimensional time 
dependent heat conduction problems: for the temperature transients in 

an infinitely long square bar; for the temperature transients in an 
infinite medium with a circular hole; and for the transient temperature 
response of a hollow cylinder within an infinitely thick wall. All 
three of these exact solutions were used as a means of establishing the 
accuracy for the numerical solution of the time dependent conduction 
equation • 

The exact solution developed by Jakob for the infinitely long 
square bar was for the case where the outer surface was suddenly 
subjected to a step change in temperature. Table 3-3 compares the 
numerical results, in dimensionless form, with those presented by Jakob. 
For the numerical solution, a diffusivity of 0.929 square meters per 
hour (10 square feet per hour) was used with a 0.2286 meter square bar. 
The temperature field was initially at a uniform value of -12.222 
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(Jo) 3aniVU3dW31 


Fig. 3-4: Code Verification Results - Multiple Zones Steady State 
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Table 3-3: Code Verification Results - Center 

Temperature Response of an Infinitely 
Long Square Bar 


at 

s 

fl o 1 

e i 

Jakob (23) 

6 i 

Numerical 

% error 

0.032 

1.000 

0.998 

0.1 26 

0.080 

0.951 

0.935 

1.600 

0.1 00 

0.901 

0.893 

0.849 

0.1 60 

0.715 

0.71 1 

0.530 
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degrees C (10 degrees F). Three minutes and fifty seconds of CPU time 
on a NAS 6650 was required to produce the results on a 17 x 17 grid, 
with convergence on space of 0.00001 and on temperature of 0.001. 
Convergence was defined as the absolute value of a simple difference of 
the old and new iterated values* 

The exact solution for the case where there was a step change in 
tenperature at the surface of a circular hole in an infinite medium was 
also solved by Jakob. Obviously, it is not possible for the present 
computer code to simulate an infinite medium. To approximate this case, 
concentric cylinders were used with an outer radius of 3.048 meters 
(10.0 feet) and an inner radius of 0.3048 meters (1.0 foot). The 
numerical solution will be valid up to the point where the initial 
thermal wave reaches the outer boundary. For the numerical problem, a 
diffusivity of 0.9290 square meters per hour (10.0 square feet per hour) 
was used. The field was initially set at 37.78 degrees C (100 degrees 
F), and the surface of the hole was suddenly raised and maintained at 
537.78 degrees C (1000 degrees F). The results obtained are displayed 
in Table 3-4. One minute and 57 seconds of CPU time on a NAS 6650 was 
required using an 1 1 x 1 1 grid with tolerance on space of 0.001 and on 
temperature of 0.01. Convergence was defined as the absolute value of a 
simple difference of the old and new iterated values. Note that even 
with the relatively coarse mesh and loose tolerances (as compared to the 
previous problem) for this single region problem, there is excellent 
agreement with the analytical results* 

The exact solution for the transient temperature response for a 
hollow cylinder within an infinitely thick wall was solved for the case 
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Table 3-4: Code Verification Results - Heat 

Flux at Surface of a Circular 
Hole in an Infinite Medium 


at 

i 

s 

q s 

kL6 s 

Jakob (23) 

q s 

kL 6 S 

Numerical 

% error 

0.01 

38.51 

38.79 

-0.72 

0.1 

1 4.1 3 

14.53 

-2.83 

0.6 

7.29 

7.40 

-1.53 

1 .0 

6.1 8 

6.28 

-1 .61 

2.0 

5.03 

5.1 3 

-2.00 

3.0 

4.49 

4.61 

-2.00 
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with the cylinder having a 1.0 foot radius; • To approximate the infinite 
plate, two concentric cylinders were used with the outer radius set at 
10.0 feet. The numerical solution should be unaffected by this false 
outer boundary up to the point where the Initial thermal wave reaches 
the outer boundary. The thermal diffusivity used was 10.0 square feet 
per hour. The temperature of the plate was initially set at 100 degrees 
Farenheit, with the temperature of the inner radius wall suddenly raised 
to 1000.0 degrees Farenheit at time equal to zero. Figure 3-5 compares 
the predicted temperature just inside the hole to the exact value 
presented in Jakob. The computed values are for an 11 x 11 grid with a 
time step of 0.001 hours. Approximately two minutes of computing time 
were needed to advance the problem 1000 time steps on an IBM 4341. 
Considering the coarsness of the mesh usee, the results are 
outstanding. 

D. Two Zones, Time Dependent 

Jaeger [24] developed an analytical solution for the time dependent 
temperature distribution in a two-layered circular cylinder. The two 
concentric layers were in perfect contact and had different thermal 
properties. Initially the temperature throughout both layers is 
uniform. At some instance in time, the inner surface of the inner layer 
is subjected to a step change in tenperature. The temperature of the 
external surface of the outer layer is held at the initial temperature. 
Jaeger did not consider internal heat generation in the solution. 

The closed form solution presented by Jaeger is somewhat tedious to 
apply* It involves a series solution that necessitates the use of 
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eigenvalues. The eigenvalues needed are the roots of a somewhat complex 
algebraic statement containing Bessel functions. Berger [25] developed 
a computer program that numerically implements the closed form solution 
developed by Jaeger. Berger's program was run for a case solved by the 
transformed time dependent conduction equation. This comparison is thus 
a check for the time dependent algorithm aid boundary condition 
algorithm for two thermally dissimilar layers. 

For the problem solved, the inner circular layer had a conductivity 
of 0.2014 x 1 0 ~ i BTU per inch-second-degrea F, a diffusivity of 0.4 
square inches per second, an inner radius of 2.0 inches and an outer 
radius of 5.0 inches. The outer layer had a conductivity of 
0.3278 x 10 -4 BTU per inch-second-degree F, a diffusivity of 0.2 square 
inches per second, an inner radius of 5.0 inches and an outer radius of 
9.0 inches. The initial temperature was 100 degrees F, with the 
internal surface of the inner layer being suddenly raised and maintained 
at 1000 degrees. These particular properties were chosen to duplicate a 
problem solved by Berger. 

Table 3-5 shows a comparison of the temperature rise with time at 
the interface of the two circular layers. Two comparisons are shown, 
one with both layers being represented by an 1 1 x 1 1 grid, and one with 
both layers being represented by a 31 x 3 grid. As can be seen, when 
both layers are represented by the finer grid, the agreement is 
excellent, with a maximum error of 3.08 percent occurring 3.5 seconds 
into the problem. The error decreases with additional time into the 
problem because the rapid temperature inc eases begin to slow down, as 
the problem advances toward the steady state solution. It should be 
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Table 3-5: Comparison Between Exact and Numerical Solutions 

of the Interface Temperature of a Circular Cylinder 
with Two Layers 


Time 
(sec ) 


Interface Temperature 
( °F ) 


Exact 

11 X 11 
Grid 
in each 
layer 

% 

Error 

31 x 31 
Grid 
in each 
layer 

% 

Error 

0.0 

1 00.00 

1 00.00 

0.00 

1 00.00 

0.00 

0.5 

1 00.00 

1 00. 1 1 

0.1 1 

100.01 

0.01 

1.0 

100.75 

1 03.06 

2.29 

101.05 

0.30 

1.5 

105.91 

113.22 

6.90 

1 07.1 8 

1.20 

2.0 

1 17.1 1 

1 30.1 1 

11.10 

1 19.63 

2.1 5 

2.5 

133.03 

151.27 

1 3.71 

136.70 

2.76 

3.0 

151 .87 

1 74.57 

14.95 

1 56.47 

3.03 

3.5 

1 72.24 

198.62 

15.32 

1 77.54 

3.08 

4.0 

193.22 

222.62 

1 5.22 

199.03 

3.01 

4.5 

214.25 

246.11 

14.87 

220.39 

2.87 

5.0 

234.96 

268.83 

14.42 

241.34 

2.72 

6.0 

274.75 

311.61 

13.42 

281.33 

2.39 

7.0 

31 1.88 

350.72 

12.45 

318.44 

2.1 0 

8.0 

346.25 

386.41 

11.60 

352.67 

1.85 
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noted that although a coarse 11x11 grid w« s adequate to determine the 
transient response of a single layered circular cylinder, the grid is 
too coarse to provide accurate results with the addition of a second 
layer; a considerably finer grid is required. This problem can become 
worse as additional layers are added, requiring finer and finer meshes 
to provide acceptable accuracies. It should also be noted that the 
interface in this problem is near the center of the two fixed 
temperature boundaries. Thus the errors at this location represent the 
maximum errors in the solution fields at a given instant in time. 

E. Multiple Zones, Time Dependent 

No closed form solutions have been found in the open literature for 
treating a time-dependent problem of a body with three or more layers 
which also contains a heat source. However, numerical solutions to such 
problems exist and these can be used for code verification. 

Stallabrass [2], Baliga [4], Marano [>], and Chao [7], using 
different numerical approaches, have each solved a one-dimensional 
time-dependent problem of a six layered slab with one of the layers 
generating heat. Table 3-6 lists the dimensions, materials and 
properties of the six layers used to verify the present code. Figure 
3-6 presents the results of computations for a time step of 0.001 
seconds with a grid having 14 nodes in the 75S-T6 aluminum, 4 nodes in 
the heater, 7 nodes in the upper insulation, 9 nodes in the stainless 
steel and 33 nodes in the ice. It can be seen that the predictions are 

in good agreement with results reported earlier. 

It should be noted that to model this series of stacked slabs, the 

present code used a set of concentric circular layers with the inner 


able 3 6: Materials and Properties Used in The Multiple 

Zone, Time Dependent, Verification Problem 


Layer Material 


Thickness 
( mm ) 


Thermal 
Dif f i^sivity 
m /hr 


7 5S-T6 
Aluminum 

2.210 

0. 1 5329 

Epoxy/glass 

Insulation 

1.270 

0.00081 

Ni chrome 
heater 

0.01 2 

0.01 282 

Epoxy/glass 

Insulation 

0.254 

0.00081 

304 Stainless 
Steel 

0. 305 

0.01 394 

Ice 

6.350 

0.0041 3 
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layer having a large radius. Some of the small discrepancies in the 
computed results are believed to be due to these geometrical effects. 

The computational time needed to obtain convergent results on a NAS 6650 
varied from two to ten minutes depending on the particular problem 
considered. To obtain the results shown, minimum tolerances of 0.00001 
on space and 0.0001 on temperature were required. Convergence was 
defined as the absolute value of a simple difference of the old and new 
iterated values. 
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CHAPTER 4 TRANSFORMATION OF PHASE CHANGE PROBLEM 

Overview of Phase Change Treatments 

Solutions to phase change problems have long been a topic of study 
in the technical community. The characteristic feature of any phase 
change problem is the coupling of two temperature fields with a moving 
boundary that not only separates the two fields, but propagates through 
them. The propagating phase front makes the problem non-linear. Latent 
heat effects and changes in thermal conductivity between phases increase 
the non-linearity. Because of the non-linearity problem, only a few 
exact solutions have been developed, and there are for very restrictive 
conditions . 

Neumann's exact solution to this problem, presented in Carslaw and 
Jaeger [26], was the first known successful solution of a phase change 
problem. Neumann solved the problem for one dimensional phase change in 
a semi-infinite region. The material was in tially at or above the 
fusion temperature, and suddenly experiences a step decrease in 
temperature at the boundary* 

Lin [27] has modified Neumann's development to obtain an exact 
solution for a quasi one-dimensional problem in a region with a varying 
cross-sectional area. To incorporate varying areas into the analysis, 
Lin simplified Neumann's approach by assuming that the entire region 
was at the fusion temperature, and by developing the relevant equations 


to solve for the interface velocity as a 


function of time. These 
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simplifications readily permitted the insertion of a transformed 
position coordinate functio into the governing equations. 

Cho and Sunderland [28] have extended Neumann's exact solution to 
allow for variable thermal conductivity in both phases. It was assumed 
that the conductivity varied linearly with temperature. The effect of 
the phase change speed with conductivity variation was investigated in 
some detail. 

Solution of the classical phase change problem, as described above, 
involves the calculation of the phase front for specified boundary 
conditions. This problem may also be solved by prescribing the location 
of the phase front as a function of time, and then determining the 
boundary conditions needed as a function of time to produce the 
specified phase front movement. This is known as the inverse problem. 
Rubinsky and Shitzer [29] have developed the exact solution to the 
one-dimensional inverse problem for any arbitrary function describing 
the movement of the phase front. Previous solutions to the inverse 
problem were for a specific function. 

Gutman [30] has constructed an approximate analytical 
one-dimensional solution that attempts to account for the effect of 
superheating or supercooling on the movement of the phase change front. 
The resulting analytical form was found by matching the inner and outer 
solutions to the problem. His method provides accurate results for a 
very limited range of relevant dimensionless parameters. 
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To date, there have been no exact solutions developed for phase 
change problems where the phase change front is a function of two space 
coordinates . 

Effort has been devoted recently towards the obtaining of 
analytical solutions by approximate methods • Ku and Chan [31] have 
developed an artificial initial condition that accounts for the 
temperature profile discontinuity at the phase change front. The 
artificial condition permits solutions to she resulting phase change 
equations by using inverse Laplace transforms. Depending upon which 
side of the phase front line the artificial initial condition is 
applied, a temperature distribution solution can only be obtained for 
one of the two phases. The technique has not been used to solve a 
two-dimensional problem as yet, but the analytical results applied to 
one-dimensional problems compare well with existing solutions. 

Zhang, Weinbaum and Jiji [32] have developed an approximate 
analytical solution to a three-dimensional time dependant phase change 
problem. The solutions are limited to very small Stefan numbers. 

Results are presented for a buried pipe example. They combine a 
quasi-steady approximation with a virtual :r ree surface method to obtain 
an axisymmetric solution for the region around the pipe wall. 
Singularities in the differential equations along the pipe would not 
permit solutions close to the pipe wall. 

The analytical and approximate analytical solutions described 
above result from the solution of partial differential equations. 

There is another class of approximate techniques that are based on an 
integral formulation of a heat balance across the melt line. One of the 
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earliest published works using this approach was by Goodman [33]. 

Goodman developed a number of approximate solutions for the 
one-dimensional problem for typical boundary conditions. Results were 
found to compare well with existing solutions when the assumed 
temperature profile was a quadratic function. 

More recently, Wang and Perry [34] have applied this technique to a 
one-dimensional problem with initial superheat. This problem is 
somewhat more involved, as there are two interfaces, one for the phase 
change front and one for the superheat line. Excellent agreement was 
obtained between the approximate solutions and a one-dimensional finite 
element code constructed to duplicate the problem studied. 

Virtually all of the available analytical models, whether exact or 
approximate, are sufficiently complex so that in most cases an alternate 
procedure must be adopted. This, coupled with the fact that most of the 
analytical models have demonstrated accuracy only for one-dimensional 
problems, strongly indicates that any feasible two-dimensional solution 
must be sought by alternative means. 

Because of the nature of the phase change problem, the obvious 
approach to take is a numerical implementation. By this method, 
generally a fixed, uniform grid is laid over the problem domain. Three 
equations are written. Two of these are conduction equations, one for 
each phase. The third equation is an energy balance, incorporating any 
latent heat effect, written across the phase change front. This 
equation ties together the two conduction equations. The advantage to 
this formulation is that the conduction equations have smooth, 
continuous derivatives of all primative variables in their respective 
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domains. This obviously enhances the likelihood of a stable, 
convergent numerical scheme. The primary disadvantage is in the 
interpolation required at the phase front. As the problem advances, for 
each increment of time, the phase front must be relocated. In most 
instances the front will be between node po nts. The diffusion 
algorithms near the phase front then need to incorporate derivatives 
using the appropriate "shortened" lengths. Because of the numerous 
interpolation calculations required, the computer codes are somewhat 
slow. Most of the current work using this approach is in the direction 
of developing faster diffusion algorithms, or making assumptions about 
the thermophysics of the problem in order to reduce the number of 
interpolati ve calculations. Lazaridis [35] for example, obtained 
two-dimensional solutions by treating the motion of the fusion front as 
quasi-one-d imensional. Thus, the temperature gradients need be taken in 
one direction only. This substantially red ices the number of 
interpolations needed. What is more, the results show satisfactory 
agreement with existing solutions. 

More recently, finite element methods lave been used to solve the 
moving phase front problem on a fixed grid. O'Neill [36] developed an 
algorithm that can be used with standard finite element heat conduction 
codes, using linear interpolation to locate the phase front within 
elements. Yu and Rubinsky [37] have treated the two relevant conduction 
equations and the interfacial energy equation as independent governing 
equations on a finite element mesh. The results of their computer 
solutions agree well with other closed form solutions. 

But even the finite element codes tend to be somewhat punitive with 
respect to code execution time, because of the interpolation/iterations 



needed to track the moving melt front. Thus, numerical techniques have 
been developed that have attenpted to either circumvent, or completely 
eliminate, the need for the interpolation caused by the phase front. 

Numerical researchers have developed a number of ways of 
circumventing the need for interpolation. These methods essentially 
involve adapting the grid to the moving melt front. Prusa and Yao [38] 
have solved a two-dimensional phase change problem of melting 
around a horizontal cylinder. The problem is solved numerically in 
cylindrical coordinates, using a radial spacial transform in the 
governing equations that essentially "stretches” or "contracts" the grid 
on either side of the phase front. The diffusion equations are then 
applied to the "stretched" regions. Duda, et al. , [39] present a 
similar technique, calling it boundary immobilization. They also use a 
stretching transformation, but in their analysis the location of the 
phase front in transformed coordinates remains at a fixed location. 

The effect, however, is essentially the same. The regions on either 
side of the phase front are "stretched" or "contracted" to generate the 
prescribed regions in transformed coordinates. 

Rieger, Projahn and Beer [18] custom fit a grid by use of a body 
fitted coordinate transformation. As the phase front moves with each 
time step, a new grid is numerically generated using the phase front as 
the upper boundary of one region and the lower boundary of the second 
region. They also solve the problem of melting around a heated 
horizontal cylinder. 

Lynch [40] uses essentially the same strategy as Reiger, Projahn 
and Beer [18], but applies a finite element mesh generator instead of a 
body fitted coordinate transform. 
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Although grid adaption eliminates the need for interpolation, 
additional calculations are required to generate the new grid (or 
function for the "stretching" transformation). The latter tends to 
offset the gains made by the former. Numerical techniques that 
eliminate the need to solve for the location of the phase front 
eliminate the need for not only interpolation, but also calculations for 
adapting the grid, since a rigid, uniform mesh can be used. In order to 
eliminate the need for any calculations foi the location of the phase 
front, the equations that govern the syster need to be reformulated. 

The reformulation, depending upon the strategy used, may result in the 
insertion of an additional approximation into the governing equations. 

Kikuchi and Ichikawa [41] introduce a special integral 
transformation that accounts for the discontinuity of the temperature 
gradient caused by the phase change front, even though the location of 
the front is unknown. This special integral transformation is known as 
the freezing index. The equations used foe iterating through the 
temperature fields are developed by applying the freezing index integral 
to the original set of governing equations. The algebraic set of 
equations needed are obtained by discretizing the associated variational 
form of the freezing index. Blanchard and Fremond [42] use a strategy 
similar to Kikuchi and Ichikawa, except that they use the homographic 
approximation to model the temperature discontinuity caused by the phase 
front. In doing this they solve for a variational equality, instead of 
the variational inequality solved by Kikuchi and Ichikawa. The primary 
advantage to the homographic insertion is that the originally 
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non-dif f erentiable change of temperature across the phase front becomes 
differentiable. 

In the heat conduction equation, temperature is usually considered 
the dependent variable, with time and space coordinates the independent 
variables. By an appropriate transformation of the conduction 
equations, however, one of the space coordinates can become the 
dependent variable, with temperature, time, and any remaining space 
coordinates the independent variables. By specifying values of time, 
temperature and one space coordinate (in a two-dimensional problem), the 
values of the independent space coordinate can be iterated for. By 
always specifying the same value of temperature at each time step, the 
movement, or migration, of an isotherm with time can be tracked. Thus, 
this technique is called the isotherm migration method. Crank and Gupta 
[43] were the first to apply this technique to a phase change problem in 
two dimensions. The advantages of this technique, in tracking a moving 
isotherm at the fusion temperature, are obvious. Results predicted by 
this approach have been found to be satisfactory, provided two numerical 
idiosyncracies caused by the technique are circumvented. First, an 
approximate method is used to calculate the position of the phase front 
for a short initial time into the problem. After the initial time 
interval, with an established isotherm, the isotherm migration algorithm 
can be used. Secondly, depending upon the geometry and/or the 
thermophysics of the problem involved, the independent space coordinate 
can become a multi-valued solution to the algorithm containing the 
dependent variables, where it is inappropriate. 


77 


Saitoh [44] extended the two-dimensional isotherm migration method 
from regular two-dimensional geometries to arbitrarily shaped, doubly 
connected two-dimensional geometries. Thi;; was accomplished by applying 
the same radial "s tretching" transformation as used by Prusa and Vao 
[38] to the governing equations before exchanging the dependent and 
independent variables. The numerical resu. ts obtained by this method 
showed excellent agreement with experimental data. 

Ozisik [45] discusses the use of a moving heat source in the 
conduction equations to account for the latent heat effect. He presents 
a mathematical development that explicitly casts the moving boundary 
problem into a standard heat conduction problem with a moving heat 
source. By doing this, the solution can inmediatedly be written in 
terms of Green’s function, and numerically implemented. 

Another technique for eliminating the necessity of locating the 
phase front, is use of the so-called high heat capacity method. In this 
method, the latent heat effect of the phase change is approximated by a 
large heat capacity over a small temperature range. This technique has 
been in existence for a number of years, and has recently undergone 
further development into more sophisticatec algorithms and applications. 
Bonacina, et al. [46] have developed a three-time level implicit scheme 
using the approach. The temperature depencent properties are evaluated 
only at the intermediate time step, thus simplifying the solution to the 
algorithms. The scheme was found to be unconditionally stable and 
convergent. 

Comini, et al. [47] present a similar analysis to the one above, 
but for a finite element code using triangular elements, so that 
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irregular geometries could be accomodated. Morgan, Lewis and 
Zienkiewicz [48] have improved this code by applying the technique with 
quadratic isoparametric finite elements. 

Hsiao [49] modified the manner in which the specific heat was 
calculated. He accounts for the latent heat, and determined the 
physical conditions of the node, by using a linear interpolation of the 
surrounding nodal temperatures. Using this linear interpolation, 
excellent results are obtained for both one and two-dimensional phase 
change problems with relatively large time steps and coarse mesh. 

Uchikawa and Takeda [50] have applied the high specific heat method 
to the irregular geometry of a casting mold. In this analysis a 
transform is used to turn the irregular regions into regular, evenly 
spaced computational zones through the use of a body fitted coordinate 
procedure. The governing equations, with the high specific heat method, 
are then transformed and solved in the transformed plane. 

One of the more recent techniques for eliminating the need to 
calculate the location of the phase front is the enthalpy method. In 
this method, the specific heat is combined with temperature to form an 
enthalpy variable in the time dependent term in the diffusion equation. 
Thus, both temperature and enthalpy are dependent variables. All of the 
earlier two-dimensional models using this technique assumed that the 
curve relating enthalpy to temperature had a finite slope. This is only 
true, however, for materials that undergo a change of phase over a 
temperature range. Shamsundar and Sparrow [51] have developed a 
variation to the earlier two-dimensional models. Depending upon both 
the temperature and the enthalpy of a given node, either temperature or 
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enthalpy is the dependent variable. When the condition of a node is 
such that a single phase exists, temperature is the dependent variable. 
When the enthalpy of a given node is be twee . the values of either phase 
(a change of phase is occurring), the tempe rarure is known, making 
enthalpy the dependent variable. This approach easily permits solutions 
for substances whose change of phase occurs entirely at a single 
temperature, as well as for those whose change of phase occurs over a 
range of temperatures. Shamsundar and Sparrow use an implicit finite 
difference scheme to numerically implement the governing equations. 

They apply the technique to PCM, a type of wax, that has relatively 
little change in thermal conductivities between phases. 

Crowley [52] applies an explicit finite difference scheme to an 
enthalpy formulation that is equivalent to that of Shamsundar and 
Sparrow. Crowley applies his algorithm to the solution of Saitoh's 
problem [44], with water as the phase change substance. Even with the 
large difference in thermal conductivities that water exhibits, the 


numerical results agreed well with published experimental data for the 
problem. 

One of the drawbacks of the enthalpy formulation, particularly for 
the case where the change of phase occurs at a single temperature, is 
that a plot of temperature against time fcr a given node tends to 
exhibit a “plateauing" effect or tendency. This arises out of the fact 
that a computational grid models or represents a discrete region in 
space. Obviously, it requires a finite amount of time to melt a 
discrete region. As a consequence of a node being held fixed at a 
single temperature for a discrete amount of time, the effect is also 
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felt by the surrounding nodes, causing a plateauing of the temperature. 
Vo Her and Cross [53] developed a smoothing technique that can be 
applied to a final set of numerical results. This smoothing technique 
has the effect of bringing the time-temperature history of a given node 
into excellent agreement with other published results. 

Schneider and Raw [54] have developed a modified enthalpy model 
that is capable of efficiently solving problems where multiple phase 
change interfaces exist. For these types of problems their modification 
reduces computational times by an order of magnitude. Thus, the 
implicit scheme used here is only applied to one-dimensional problems, 
but they indicate that it can be easily extended into two dimensions. 

Voller [55] provides an alternate method of discretizing the 
enthalpy formulation. In this approach, the sensible and latent 
heat terms are discretized separately, thus isolating the non-linearity 
of the problem as a nodal latent heat source term. His implicit 

finite-difference scheme yields computational savings of twenty to fifty 
percent • 

Tacke [56] has developed a formulation of the enthalpy method which 
removes the "plateauing" effect in the time-temperature history of a 
node. This was accomplished by applying linearized temperature profiles 
near the phase front for those nodal control volumes containing the 
front. His numerical solutions substantially reduced the plateauing, 
while showing excellent agreement with analytical solutions. The 

linearization does increase the computational time required, but only 
slightly. 
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Recently, Schneider [57] has used the enthalpy form of the energy 
equation, and, in the melted liquid region, coupled it with the momentum 
equation to account for the effect of free convection effects on the 
movement of the phase front line. The prollem is solved in a 
rectangular region, with a rise in temperai are occurring at only one 
vertical wall. Boucheron and Smith [58] have solved essentially the 
same problem, but they couple the momentum equation throughout the field 
for both phases. They specify for the sol:d phase a very high value of 
viscosity in order to insure that the veloc ities arising from the 
solution of the momentum equations in this region are negligible or 
zero. This approach allows a more mathematically straightforward 
treatment of the phase front. 

Because of the relative ease of formulation, simplicity in 
numerically discretizing the resulting equations, and proven accuracy 
and stability, the enthalpy formulation was selected to solve the phase 
change portion of the current problem. The irregular region is mapped 
into a rectangle using a body fitted coordinate spacial transform. The 
governing equations must be similarly transformed, and solved on the 
transformed grid. Writing a general "numerical" energy balance in the 
transformed plane is virtually impossible. As was discussed in the 
previous chapter, the square, regular computational cells in the 
conputational or transformed plane unmap into nonuniform, non-orthogonal 
cells of varying area in the physical plant. Consequently, the 
transformed partial differential field equations were differenced rather 
than developed by conservative principles. Because thermal conductivity 
is not a constant for this problem, there are two possible forms of the 
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field equation , the conservative form and the non - conserva tive form, 
which are mathematically equivalent. When treated numerically in the 
same fashion, the result will be two algorithms that are not equivalent. 
The rationale for selecting the equations type is the topic of the next 
section. 

Importance of Equations Type - Neumann Comparison 

Since the energy balance at the phase change front is essentially 
"buried" in the weak formulation of the phase change equations, the 
discretized form of these equations can become important for an accurate 
solution. One important distinction in equation form is conservative 
versus non-conservative. The term "conservative" means that the 
"purest" form of an equation is preserved or "conserved"; while 
^ ons e rva ti ve means that some change has been made to the "pure" form 
(i.e#, differential operators have been carried through, etc.). This 
difference in form for the enthalpy formulation is shown in detail 
below. 


Figure 4-1 depicts the two zones that need to be considered in the 
traditional formulation of a phase change problem. The equations that 
apply $ H only thermal conduction with no density change is considered, 
are: 


in liquid 
region 



(4-1 ) 


in solid 
region 


(4-2) 



83 



Figure 4-1 The Two Zones needed in the 
Traditional Formulation of 
the Phase Change Problem 
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k s k L — = P s v n A (4-3) 

3n 3n 

on f (x, y , t) =0 

Tj_(x,y # t) = T f i = s and L (4-4) 


where n - outward normal to the interface, into 

the liquid 

v n = velocity of the interface in the normal 
direction 

f(x,y,t) = function describing the interface 

separating the solid and liquid regions 

The above two field equations, coupled by an interfacial energy 
balance, Eq. (4-3), can be reduced to a single non-linear field equation 
that eliminates the need for any computation of the location of the 
phase front. This is known as the weak formulation, and for phase 
change problems is called the enthalpy formulation. The term "weak 
formulation" means that less direct information is provided in the 
solution of the equations (in this case, the location of the phase 
change line). Thus, the weak formulation for this problem becomes: 


9t\ 3 / 3t\ 

— + — ( k — J (4-5) 

9x/ 3y \ 9y / 

H 

when H < H sm , T = (4-6) 

and k = k s ; 



when H sm < H < H Lm , T = T m 

H-H sm 

and k = k s + (k L - k s ) 

H Lm -H sm 


(4-7) 
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H Lm- H 


and when H > H Lm , T = T n 


(4-8) 


PC L 


and k = k Ti . 


The above set of equations is much more easily treated computationally 
than the previous set. The location of the melt line can be determined 
indirectly from the distribution of enthalpies on the solution domain. 

Equation (4-5) is written in the conservative form. If the 
operators are carried through, Eqs. (4-5) becomes 


9h 

3 T 

3^T 

9k 

9T 

9k 

9t 


— = k 

+ k 
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-- + 
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(4-9) 

9t 

3x^ 

3/ 

9x 

9x 

dy 

dy 



and is now cast in a non-conservative form. 

Written for the one-dimensional case, Eqs. (4-5) and (4-9) become, 

respectively 
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/ 
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(4-10) 

3t 

9x 
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3H 
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_ _ 

= k 

_ , 

+ -- 

— _ 

(4-1 1 ) 

9t 


3x* 

9x 

9x 



If a simple explicit differencing procedure is applied per the 
general one-dimensional grid depicted by Fig. (4-2), Eqs. 4-10 and 4-11 
become 
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H k+1 _ H k 


At 
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-- 2 

(Ax) 


k+1 

k j + 1 + 
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k+1 

kj k+1 
— (Tj + 1 


k+1 


k+1 k+1 

kj + k j-1 


k+1 

(T j 


k+1 

h-i 


(4-12) 


2 
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H k+1 _ H k \ 1 

At J ( Ax ) 

k+1 k+1 k+1 k+1 . (4-13) 

( k j + 1 ” k j-l)(Tj + 1 “ \ 

+ ““ ) 

where the superscript k denotes the time step and the subscript j the 
nodal location. The algorithms obtained by using Eqs • (4-12) and (4-13) 
are clearly not equal to each other since thermal conducting varies as 
the phase changes. Error is introduced into any mathematical system 
when that system is numerically discretized. The errors introduced by 
the algebraic expressions in Eqs. (4-12) and (4-13) are obviously 
different, with one contributing potentially more error than the other. 

As the difference between conductivity increases between phases, so 
does the potential disparity between Eqs. 4-12 and 4-13, especially for 
a computation with a node undergoing a charge of phase. For this case 
the melting node would be bounded by at least one all liquid node and at 
least one all solid node. It should be noted that for constant 
conductivity between phases, Eqs. (4-12) and (4-13) reduce to 
identically the same algorithm. The problem arises only with a change 
in conductivity with phase in the discretized equations. 

It should also be mentioned that in the traditional formulation for 
phase change, this problem cannot occur, regardless of whether the 
conservative or non-conservative form of the equations are differenced, 
since conductivity is constant in a given zone. The problem arises in 
the weak formulation of the equations where the location of the 
discontinuity is "buried” within the field equation. 



k+1 k+1 k* 1 k+1 

kj (T j + 1 - 2Tj + Tj.-,) 
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Equations (4-12) and 4-13 can be checked for accuracy by comparison 
with an approximate analytical solution developed by Neumann and 
presented in Lunardini [59], 

Figure 4-3 depicts the problem solved by Neumann. The problem can be 
formulated as: 


1 9T S 

9x 2 3t 


(4-14) 


1 9T L 

9x * Ol a t 


(4-15) 


Limit ( Tl ) = (4-16) 

x + 00 

Tg (0,t) = T sur jp ace (4-17) 

where T SU rface > T f 

at the interface 

T s (x,t) = T l (x,t) = T f (4-18) 

and 


9 T S 


9x 


— 

- k L - 

pX -- 

(4-19) 



3t 



In the solution of Eqs. (4-14) through (4-19), if the initial 
tenperature of the liquid is Tf , and the melt distance into the problem 
is small, the solution can be approximated by: 


■-V- 


(4-20) 


2k s (T f -T surface )t 
Tables 4-1 and 4-2 show in detail the numerical computations for 


solving the one-dimensional phase change problem using the explicit 


Solid - Region 1 




the Phase Change Problem 
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EXPLICIT ALGORITHM FOR 1-D CONSERVATIVE EQUATION: 



0.0 

0.00 

0.00 

0.00 

1.4160 

1.4160 

0.1 

0.00 

26.04 

26.04 

1.4160 

1.2164 

0.2 

26.04 

23.05 

49.09 

1.2164 

1.0398 

0.3 

49.09 

20.39 

69.48 

1.0398 

0.8834 

0.4 

69.48 

18.05 

87.53 

0.8834 

0.7450 

0.5 

87.53 

15.98 

103.51 

0.7450 

0.6227 

0.6 

103.51 

14.14 

117.65 

0.6227 

0.5143 

0.7 

117.65 

12.51 

130.16 

0.5143 

0.4184 

0.8 

130.16 

11.08 

141.24 

0.4184 

0.3335 

0.9 

141.24 

9.80 

151.04 

0.3335 

0.3200 


Melt Time for this Node Per 
the Neumann Solution: 0.7456 sec. 


Error (using interpolated time): 7.37‘> 


Table 4-1 Numerical Computation for the Conservatively 
Differenced Phase Change Equation 
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EXPLICIT ALGORITHM FOR 1-D NON-CONSERVATIVE EQUATION 


,k+l 


yk i At rp k (T k -?T k +T k i 

= H + 7UZY Lk j ( Vi j 


+ 


“4 


] 


Time 

H k 

AH H k+1 

K k K k+1 

(sec) 

(BTU/lbm) 

(BTU/hr-ft-°F) 


0.0 

0.00 

0.1 

0.00 

0.2 

34.26 

0.3 

60.64 

0.4 

80.95 

0.5 

96.59 

0.6 

108.64 

0.7. 

117.92 

0.8 

125.07 

0.9 

130.57 

1.0 

134.80 

1.1 

138.06 

1.2 

140.58 

1.3 

142.50 


0.00 

0.00 

34.26 

34.26 

26.38 

60.64 

20.31 

80.95 

15.64 

96.69 

12.05 

108.64 

9.28 

117.92 

7.15 

125.07 

5.50 

130.57 

4.23 

134.80 

3.26 

138.06 

2.52 

140.58 

1.92 

142.50 

1.49 

143.99 


1.4160 

1.4160 

1.4160 

1.1534 

1.4160 

0.9512 

0.9512 

0.7955 

0.7955 

0.6567 

0.6567 

0.5833 

0.5833 

0.5122 

0.5122 

0.4574 

0.4574 

0.4153 

0.4153 

0.3828 

0.3828 

0.3578 

0.3578 

0.3386 

0.3586 

0.3238 

0.3238 

0.3200 


Melt Time for this Node per 
the Neumann Solution: C.7456 sec. 


Error (using interpolated time): 74.35% 


Figure 4-2 


Numerical Computations for the Non-Conservati vely 
Differenced Phase Change Equation 


L 
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conservative and non-conservative algorithms. The calculations are 
performed for water as the phase change substance, and are done to 
determine the required melt time for the first node only. Figure 4-4 
provides a schematic of the problem solved. Tables 4-1 and 4-2 provide 
further insight into the nature of the numerical error that is 
introduced into the solution. 

For the algorithm based on the conservative form, the heat flux 
must be written at each edge of the node. To estimate the thermal 
conductivity at the node edge, the average is taken of the 
conductivities of the two nodes having a common edge. The conductivity 
of the node j is determined by a linear interpolation between the 
solid and liquid conductivities, depending on the percentage of the node 
that has melted. For the algorithm based on the non-conservative form, 
no averaging is needed since the first order derivatives are evaluated 
based on values of the two adjacent nodes, which never change until the 

j node needs to be updated, depending on the percent of the node 
that has melted. 

As the node "warms up" and melts, both solutions track each other 
with good agreement. The solutions begin to diverge, however, during 
the melting of the last third of the node. The non-conservative 
differencing underpredicts the heat flux entering node "j". Clearly, 
approximating the conductivity derivative, or, the change of 
conductivity through three nodes, is less accurately modelled in 
non-conservative differencing. Averaging the conductivities between 
nodes for the edge, and using these averages in a flux derivative, is a 
much better numerical approximation. Note that the relative error 
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AX = 


Figure 



k , = 1.416 BTU/hr-f t-°F 

'Olid 

T = 32°F 

c. 


*Sol id — ^ — ^liquid 
T = T m 


k, - .,= 0.320 BTU/hr-ft-°F 

iquid 


node 


= T + 30°F 0 t 


0 


4-4: Schematic, with Properties, of One-Dimensional 

Grid Used to Obtain Solutions Presented in Tables 
4-1 and 4-2. 
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between the two methods, just in a first node calculation for this crude 
discretization, differ by an order of magnitude* Using a smaller mesh 
or finer time step does not improve the approximation appreciably* 

Figures 4-5 and 4-6 show the numerical results of Eqs. (4-10) and 
(4-11) for a Neumann class of problems* The solutions plot the movement 
of a melt front with time. At the start of the problem, all the nodes 
are solid and at the fusion temperature. Suddenly, one edge of the 
region experiences a step change in temperature, and is held constant at 
that temperature. The substance under consideration is water, which has 
a difference in conductivity between the two phases of approximately a 
factor of four. Solutions are plotted for step changes in the edge 
temperature of 10, 20, 30, 40 and 50 degrees Farenheit above the fusion 
temperature . 

For this problem, which solves for the melt front passing through a 
number of nodes into the grid, the Crank-Nicholson implicit 
finite-difference scheme was employed in the numerical implementation of 
Eqs. (4-10) and (4-11), This scheme was chosen for the two-dimensional 
problem because it is unconditionally stable for the conduction 
equation, [60]. The solutions for the conservative and non-conservative 
forms with Neumann's solution using Crank-Nicholson differencing 
provides further justification for choosing the correct form for the 
two-dimensional problem. Using the Crank-Nicholson differencing scheme, 
Eqs. (4-10) and (4-11) for the problem solved become: 



Melt Front Distance (inches) 









Melt Front Distance (inches) 


96 
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Figure 4-6: Comparison of Neumann's Solution to the Enthalpy 
Equation Differenced in Non-Conservative Form 
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All second order derivatives are approximated using three point 
differencing; all first order derivatives using two point differencing, 
Gauss-Seidel iteration was used so that the temperature profile with the 
effect of the phase change could be obtained at any time step. 

As can be seen from Fig. 4-5, Cr ank-Nicholson differencing of the 
conservative form of Eq . (4-5) yields results that compare favorably 
with Neumann's approximate analytical solution. However, the agreement 
deteriorates with increasing values of the step change temperature. 

This is to be expected, since if all numerical criteria for a 
computational zone remain constant, a more rapidly changing solution 
with real time is apt to be less well approximated by the derivatives. 

The large error introduced into the solution by differencing the 
non-conservative form of Eq. (5-4) is very apparent in Fig. 4-6. The 
error (defined as the difference between Neumann's value and the 
numerical value divided by Neumann's value) is fairly consistent, and is 
equal to approximately 100 percent through the entire solution. The 
numerical solution shows about twice the amount of time needed to 
achieve a desired melt distance than that predicted by Neumann's 
solution . 

In the numerical simulation of the weak formulation of the phase 
change equations, attention needs to be given to the conservative and 
non-conservative forms of these equations. For the particular problem 
at hand, clearly the conservative form should be used. 

Transformation of Phase Change Equations 

Since the phase change equations will be solved in the transformed 
plane, those portions of the equations that have a spacial dependency 
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must be transformed. Because it has been ;hown that the conservative 
form of the field equation must be used, only the transformation 
operators for the first order derivatives need to be used. From Chapter 
2, the appropriate derivative operators needed to effect the 
transformation of the field equation are: 

3 1 / 3y 9 3y 3 

3x J \3n 35 35 Bn 


(4-25) 


3 1 / 3x 3 3 

3y j V 8 5 Bn 3g 35 


Applying these to Eq. (4-5) yields: 
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Expanding Eq . (4-28) produces the following: 
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Thus Eq. (4-31) may be written as 
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Equation (4-35) is the transformed enthalpy equation [Eq. (4-27)] that 

must be solved in the transformed plane. 

The boundary condition equations are th ; same as those for a 
layered conduction problem, that is, equality of tenperature and 
heat flux at the interfaces: 

I 

T (4-36) 
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For a convection condition: 
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s 


Transformation operators are needed to provide a value of the 
temperature derivative normal to a surface at a specified point. Again 
from Chapter 2, these are 

9 1 / 9t 9t 

-- = I y 6 — 

dn jSy~ \ an 

To a line of 
constant n 



On the presumption that the boundary is along the upper or lower edge in 

the transformed plane (a line of constant n), E^s. (4-37) and (4-38) 
become 



Equations (4-41) and (4-42) are the transformed normal heat flux 
equations that must be solved in the transformed plane. 

Discr etization of the Phase Change Equations 

Figure 4-7 depicts a section of the grid in the transformed plane. 
Because of the fact that the individual elements in the real plane may 
be irregular and unsymmetric, all eight nodes surrounding the node of 
concern must be involved in the differencing. It becomes more apparent 
as to the manner in which Eg. (4-31) needs to be differenced if the 
equation is "numerically" broken down as follows: 




Figure 4-7 : 


Schematic of Grid with a General Node 
Assignment in the Transformed Plane 
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(4-43) 


As can be seen, the Crank-Ni cholson numerical scheme is being used 
Three point differencing is applied to all second order derivatives and 
two point differencing to all first order derivatives. Note that the 
locations (1+1/2, J), (1-1/2, J), (I,J + 1/2), (I,J-1/2) are at the 
boundaries of the central node. If the central node should be 
undergoing a change of phase, it will have at least one solid and at 



least one all liquid node adjacent to it* Assigning a correct value of 
conductivity at the "edge" now becomes a concern. The approach used by 
Marano [5] and Chao [7] will be used. The conductivity will be 
estimated by simply averaging the conductiv; ties of the two nodes having 
this edge in common. The differenced field equation for the phase 
change region then becomes 
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(4-44) 


Note that with the insertion of "M" and " 1-M", the algorithm can become 
anything from purely implicit (value of "M" equal to zero) to purely 
explicit (value of "M" equal to 1.0). For an "M" of 0.5, the algorithm 
becomes the Crank-Nicholson scheme. 


The boundary condition for equality of normal heat flux at a solid 


boundary becomes: 



(Ti,2 - Ti,l ) - 







( ^i , j max 

(Ti , j nax “ T i , j max- 1 ) 

, jmax , j max 

^i,jmax T i + 1,jmax “ T i-1 # jmax 

7 ■ • / y . 1 2 

J i , j max T i , j max 

(4-45) 

The algorithm is similar for a convective boundary, except that 
h (T ( I , 1 ) - T infinity) or h(T(I,J max ) - T infinity) replaces either the 
left or right hand side of Eq.(4-45), depending on whether the 
convective boundary is on the upper or lower surface of the grid. 

It should be noted here that one may wish to account for the latent 
heat of fusion for the half-cell at the boundary in the phase change 
region for large mesh sizes. Though the mathematical formulation is 
correct, it applies only to a point. When numerically modelling a 
problem, however, the numerical equivalent applies to a discrete region 
in space. In the problem being modelled, cne half of this discrete 
region at the boundary absorbs heat in the amount of the latent heat of 
fusion. For a large discrete region (large mesh size) the error 
introduced in directly applying Eq. (4-45), without modifications to 
incorporate the latent heat of fusion in tie half cell, may be 



significant. 
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Logic for the Numerical Solution of Phase Change Field Equation 

One approach for numerically implementing Eq . (4-31) under the 
conditions specified by Eqs. (4-6) through (4-8) in a manner that is 
computationally convenient, is to write three field equations; one each 
for the solid, melt, and liquid states. In doing this, only the equation 
for the melt state needs to be written with an enthalpy term, while the 
other two equations can be left in terms of the temperature. The nodal 
enthalpy is then calculated separately from the iterated temperature. 

The three field equations become 



where H = P s c s T and k = k s at the node (4-46) 



^ode H sm 

where k = k s + (k L - k s ) (4-47) 

H Lm “ H sm 


and T = T m at the node 



where H = H Lm + P^c L (T-T m ) 


and k 


at the node 


(4-48) 



1 09 


It should be noted that the differencing for the solid and liquid field 
equations must be almost identical to the melt equation (which resulted 
in Eq. (4-44)). The only difference is that the specific heat is not 
combined with temperature in the time derivative to give the enthalpy 
variable. The reason for this is that even though a node may be all 
solid or liquid, it may have next to it a node undergoing a change of 
phase. The heat flux through a cell edge adjacent to a melting node will 
not be properly accounted for unless there is a proper "averaging" of the 

conductivities • 

Figure 4-8 graphically depicts the iteration logic for a node in a 
phase change region. A check is first made on the previously calculated 
enthalpy of a node. Based on this check, either the solid, melt, or 
liquid algorithm is used. A new value of erthalpy is calculated, and a 
check is performed to see whether this value is within the range of 
enthalpies for which the algorithm is valued. If the check is positive, 
the subroutine moves on to the next node. If it is negative, a new 
algorithm is chosen, based on the new enthaipy value. In the event that 
none of the three algorithms calculate an enthalpy value valid within its 
range, a counter in the subroutine (which increases by a single integer 
each time an algorithm is used for a given node) will terminate the 
calculations, indicating that convergence was unachieveable for a node in 

the melt region. 

It should be noted that in the melting subroutine a convergence 
problem often occurs at those points which represent the onset and 
conclusion of a phase change. In other words, for example, at the onset 
of melting, the solid phase algorithm will calculate an enthalpy that is 
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Figure 4-8. Code Logic for the Melting Subroutine 




just Slightly above its "legal" range, and the melt algorithm will 
calculate an enthalpy that is just slightly below its "legal" range. 

This convergence problem is alleviated, with virtually no effect on 
accuracy, by numerically allowing the "legal' ranges to overlap each 
other slightly. 

Appendix F lists the key subroutine tha : implements the algorithms 
for the enthalpy formation of the phase change equations. Appendix G 
lists the key subroutine that implements the conductive boundary 
condition for a phase change region. 

Comparison with One-Dimensional Iced Airfoil 

Figure 4-9 shows a comparison between results obtained using the 
computer code developed here and Marano's [5] one-dimensional code 
developed for a composite body with phase change, which was intended for 
the purpose of modeling an iced airfoil. Beth codes determined the 
transient response of a "standard" electrothermal deicer, as defined in 
Marano's work in Table 3-6. The standard deicer was initially at a 
uniform temperature throughtout. At the initial instant of time, a 
nichrome heating element was engaged having an equivalent surface heat 
flux of 25 watts per square inch. Figure 4-9 plots the movement of the 
melt line in a layer of ice as the composite "warms up", beginning with 
initial temperatures of 10 and 20 degrees Farenheit below the melt 
teirperature. In the current study, the one-dimensional problem was 
simulated by usinp two-dimensional concentric cylinders to identify 
layer boundaries. This problem, for uniform conditions in the angular 
direction and a sufficiently large circle radius compared to layer 



Melt Distance Into Ice (Inches) 




thickness, will model 


Marano’s one-dimensional problem with negl: gible error due to geometry. 
The temperatures predicted by both codes conpare favorably. 

There is initially a slight overpredi* tion in the melt front, with a 
crossing of lines and a subsequent slight i, nderprediction for later times 
in the problem. This result would be expe( ted in light of the comparison 
with the transient response of the ice/abr<ision shield interface 
temperature presented in Chapter 3, Figure 3-6, for the pure conduction 
problem. As the initial temperature rise ft the interface slightly leads 
the results shown by Baliga [4], Marano [5 , and Chao [7], it is to be 
expected that initially the melt line slightly leads in the phase change 
problem. Also, as the temperature rise lags slightly later in the 
problem, one would expect a corresponding lag in the melt line for the 
phase change problem. 
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CHAPTER 5 EXPERIMENTAL COMPARISON WITH AN ICED AIRFOIL 

Over the years extensive research in the areas of ice accretion, 
de-icing, and the effects of icing on aero-performance have been 
investigated in a subsonic icing wind tunnel located at the NASA Lewis 
Research Facility in Cleveland, Ohio. The IRT (Icing Research Tunnel) is 
essentially the same as any closed loop subsonic wind tunnel, with two 
significant differences. First, the tunnel contains a bank of chillers 
which obtain/maintain tunnel temperatures well below freezing. Secondly, 
a spray rig exists upstream of the test section to inject a water mist 
into the airstream. Under proper conditions, ice will accrete on an 
object in the test section. The intent is to simulate the natural 
accretion of ice on a body as it would occur in flight under icing 
conditions . 

Recently, a battery of tests was undertaken in the IRT to 
investigate the performance of an electrothermal de-icer pad installed in 
a section of a UH1H helicopter rotor blade. A portion of the test 
results have been reported and analysed by Leffel, et al. [9]. The 
electrothermal de-icer used in the blade was designed and manufactured by 
the B. F. Goodrich Co. The testing was conducted in four phases: dry 

air tests, wet air tests, accretion documentation tests, and, finally, 
de-icing tests. The computer code developed in the current work was used 
to simulate the thermal response for a section of the blade for one of 


the de-icing tests 
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Figure 5-1 illustrates the layered cons ’.ruction of the blade 
section that was fabricated for the test. The assembly is similar to 
that of a standard de-icer, except that the substrate is composed of 
three separate layers: an aluminum skin whi :h is wrapped around a 

doubler that sits on the D-spar. At the leading edge, the doubler 
thickens and becomes what is known as the no?eblock. The noseblock, 
typically, is made from brass, not aluminum. Each of the layers are 
bonded together with either an epoxy glue or a film adhesive. Note that 
with the bonding materials, there are a total of thirteen layers, not 
including any accreated ice. The heater was divided into eight separate 
one inch zones along the arc. Each zone car fire independent of the 
other zones. 

As the blade was being fabricated, three layers were heavily 
instrumented with thermocouples. The thermocouple placements within the 
layers are graphically depicted in Fig. 5-2. The thermocouples were 
placed at the inner side of the D-spar, the inner side of the heater 
mat, and the outer side of the abrasion shield. The rows of 
thermocouples were placed arc-wise at the heater segment centers, 
through three "cuts" of the test section. The transient responses of 
similarly located thermocouples in the three cuts were an indicator of 
uniformity of de-icing performance in the span-wise direction. Each cut 
also serves as a backup to insure a reading at a specified location in 
the event of a thermocouple failure at a similar location. 

The material properties of the layers in the test section, along 
with average values of their thicknesses, are presented in Table 5-1. 
Modeling a de-icer pad consisting of the thirteen layers (fourteen with 
a layer of accreted ice) for a reasonable number of nodes along the arc 




Figure 5-1: Electrothermal De-Icer and Blade Construction 
of a UH1H Helicopter Rotor Blade 









Figure 5-2: Instrumentation Planes and Thermocouple 

Placements for the UH1H Blade Test Section 
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ACTUAL BLADE GEOMETRY 

NUMERICAL SIMULATION j 

Values of k 

in BTU/hr-ft 

- °F 

Values of 

k in BTU/hr- 

f t-°F 

Values of ct 

in Square ft/hr 

Values of 

ot in Square 

ft/hr 

Values of thickness in inches 

Values of 

thickness in 

inches 

Layer 

Pr operties 

Thickness 

Thickness 

Properties 

Layer 

Abrasion 

k = 8.7 

0.030 

0.030 

k = 8.7 

Abrasion 

Shield 

ot = 0.1 5 



a = 0.15 

Shield 

Adhesive 

k = 0.1 

0.0168 




Epoxy 

a = 0.0058 





Insulation 

k = 0.22 

0.0138 

0.0388 

k = 0.1 

Insulation 


a = 0.0087 



a = 0.0058 


Adhesive 

k = 0.1 

0.0082 




Epoxy 

a = 0.0058 





Heating 

k = 60.0 

0.0065 

0.0065 

o 

II 

At 

Heating 

Element 

a = 1.15 



a = 1 .1 5 

Element 

Adhesive 

k = 0.1 

0.0082 




Epoxy 

a = 0.0058 





Insulation 

k * 0.22 

0.1 38 

0.1 544 

k = 0.1 

Insulation 


a = 0.0087 



a = 0.0058 


Adhesive 

o 

• 

o 

II 

At 

0.0082 




Epoxy 

a = 0.0058 





Blade 

k = 8.7 

0.02 




Skin 

ct = 0.1 5 





Fi lm 

k = 0.1 

0.01 




Adhesive 

a = 0.0058 





Aluminum 

k = 102 

0.05 

0.265 

k = 102 

Alumi num 

Doubler 

a = 2.83 



a = 2.83 

D-Spar 

Film 

k = 0.1 

0.01 




Adhesive 

a = 0.0058 





Aluminum 

CM 

O 

t— 

II 

At 

0.1 75 




D-Spar 

a = 2.83 




l 


Table 5-1 Actual Blade Thicknesses and Material Properties vs. Those 
Used in Numerical Simulation 
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would be computationally quite difficult; the estimated CPU time would 
be extremely large. Thus, a number of layers were "lumped" together in 
the numerical simulation of the blade. The adhesive epoxy on both sides 
of the upper and lower insulation layers has been lumped into the 
insulation. The total thickness assigned : s the sum of three individual 
thicknesses, with the material properties being those of the adhesive. 
The blade skin, doubler and D-spar, with the two layers of film 
adhesive, have been lumped into a single "D-spar" layer. Again, the 
total assigned thickness to this layer is - :he sum of the five individual 
thicknesses, with the assigned properties being those of aluminum. 
Lumping the layers in this fashion for the lower insulation and the 
D-spar should have very little effect on the thermal transients in the 
heater and abrasion shield. In the blade construction the lower 
insulation has 16 times the thickness of the film adhesive, with the 
properties of the film being on the same order of magnitude as the 
insulation. Since the lower insulation is 10 times thicker than the 
upper insulation, most of the energy initially generated in the heaters 
will be driven toward the abrasion shield. Thus, any lumping below the 
lower insulation should have a negligible effect on the temperature 
transients in layers above the heater, especially for short real times 
into the problem. Lumping together the upper insulation with two epoxy 
layers, and having material values of the epoxy, may slightly retard the 
temperature at the abrasion shield. The insulation material values are 
slightly higher for both conductivity and diffusivity. 

The UH1H airfoil is the same as the NaCA 0012. Using dimensionless 
NACA 0012 coordinates provided by Abbot and Doenhoff [61], the 
coordinates for those portions of the airfoil containing the heaters 



were generated. These coordinates were assigned node numbers for 
subsequent numerical computations, which are displayed in Fig. 5-3. 

These nodes represent the coordinate locations for the outer edge of the 
abrasion shield. A coordinate generating subroutine, beginning with 
these nodes as a starting point, then generated the boundary coordinates 
for any inner or outer layers. On the outer edge of the abrasion 
shield, the nodes were spaced at one eighth inch intervals. This made 
the heater zones one inch wide, with each heater having a node with an 
adjacent heater. Nodal locations 21, 29, 37, 45, 53, 61, 69 and 77 
represent the thermocouple arc locations on the UH1H test section. 

The test case simulated by the numerical code was designated in 
Leffel's work as reading 234, position 5 (node 69), thermocouples S3 
(abrasion shield) and 26 (heater). For this particular test, the wind 
tunnel test speed was 100 mph, the ambient temperature was 16 °F , the 
angle of attack was zero degrees, the heater power density was 8 watts 
per square inch, and the air liquid water content was 2.2 grams per cubic 
meter with an average droplet diameter of 19.2 microns. The accretion 
test for this run showed that approximately three eighths of an inch of 
ice had accreted on the test section near node 69. For the first cycling 
of heater zone "G H (see Fig. 5-3), the heater was engaged for twenty 
seconds, and then disengaged for thirty seconds. 

A comparison of the thermocouple data generated at position 5 (node 
69) with the numerical simulation predictions is presented in Fig. 5-4 
for the conditions described above. As can be seen, there is excellent 
agreement between the numerical and experimental results for that portion 
of the cycle for which the heater is engaged, and for the first third of 
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KD 


Figure 5-3*. Node Numberings and Heater Designations for Numerical Model of 
the UH1H Rotor Test Section. Nodes Shown on Axes Represent 
Location of Outer Surface of Abrasion Shield 
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Figure 5-4: Comparison of an Experimental Result 
(Position 5, Reading 234) with 
Numerical Prediction 



that portion with the heater disengaged. For both the heater and 
abrasion shield, the numerical model under?) edicts the magnitude of the 
heat dissipation in the test case for the latter portion of the cycle 
when the heater is disengaged. Since the numerical simulation is 
clearly modelling the transients accurately for the "warm up" and the 
first portion of "cool down", there is obviously some physical 
phenomenon occurring in the test case that is not properly accounted for 
in the numerical simulation of the latter portion of the "cool down". 

One obvious possibility is the loss of ice (either through natural 
shedding or with an "assist" by the initial warming) at this location. 
With less ice at this location (ice is an excellent insulator) the heat 
would have dissipated much more quickly in the test case. A second 
possibility, and one that is much more likely, is that the thermal 
properties of the materials in the test section are changing with 
temperature. The numerical code can handle this condition with only 
slight modification, provided that the material properties are known 
a function of temperature. Whether or not material properties for the 
insulation, epoxy glue, and film adhesive vary as a function of 

temperature is not currently known# 

A test case of the numerical code was run to determine the extent 

of any geometric effects on temperature transients that current one-and 
two-dimensional codes cannot model. The same thickness and materials 
were used as in the comparison above, but this time the heater density 
was 16 watts per square inch. In this simulation, in which all of the 
heaters were fired, a case is provided where effects due to geometry 
only could be investigated. Table 5-2 displays the results of selected 
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Table 5-2 Time-Temperature History of Selected Interface Nodes for a Numerical 
Simulation of a UH1 H Rotor Section with a De-Icer Input of 16 W/in . 



points on layer interfaces along the leadinc edge arc (where curvature 
is the highest). As can be seen, the outer layers experience a slight 
temperature drop, while the inner layers experience a slight temperature 
rise, for those regions near the leading edqe. This is exactly what 
would be expected. As the thermal wave moves outward from the heater, 
there is more mass to absorb the generated energy. As the thermal wave 
moves inward from the heater, there is less mass, as a consequence of 
curvature, to absorb the generated energy. Consequently, one would 
expect a slight temperature drop outward from the heater, and a slight 
temperature rise inward from the heater, in contrast to results 
predicted using a one-dimensional model. For the particular case run, 
the temperature rise at the blade's stagnation point is 2.35 degrees 
Farenheit lower due to curvature/geometry. Out of a total temperature 
rise of 27.54 degrees expected for the one-dimensional case, this 
represents an error of 8.5% at this point due to geometric effects 25 
seconds into the problem. Obviously, the magnitude of the error will be 
primarily a function of distance from the source of the thermal 
disturbance, the strength of the thermal disturbance, the degree of 
curvature, and the amount of real time intc the problem, for this 


particular simulation. 
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CHAPTER 6 CONCLUSIONS 

In this thesis a computer code has been developed that is capable 
of handling any number of irregularly shaped layers in a composite body 
undergoing a transient conduction process* The code is capable 
of handling thermal generation within any of the layers, as well as 
phase change in an outer layer* The computer code was verified for a 
wide variety of test problems* 

The computer code was employed to simulate the actual transient 
thermal response of a UH1H rotor blade equipped with an electrothermal 
de-icer system* Good agreement over the majority of a heating cycle 
with test data was obtained. Additional runs with the code on the UH1H 
blade cross-section clearly showed the effect of geometry/curvature on 
the thermal transients. 

It is believed that this code can be a very useful tool in airfoil 
ice protection design. It may also be used to determine where and under 
what conditions one-dimensional codes will yield satisfactory results 
(with obvious savings in time and money), and where a two-dimensional 
code is required. 

Additional work is needed to apply more numerically sophisticated 
techniques into the spacial transform, the conduction, and phase 
change portions of the problem to accelerate convergence. The run made 
for the experimental comparison required 30 minutes of CPU time on a 
CRAY-XMP located at the NASA-Lewis Research Facility in Cleveland, 


Ohio. 



Finally, additional IRT tunnel experimentation is needed in order to 
provide additional insight into the thermc -physics and other related 
phenomena effecting electrothermal de-icer performance. 
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ALPHA VALUES FOR 

THE LAYER 

3.930280 

3.930252 

3.038534 

3.038523 

1,761444 

1 .761443 

1.038439 

1,038438 

0.621851 

0.621850 

0.377859 

0.377859 

0.232761 

0 . 232761 

0.145231 

0 . 145231 

0.091716 

0.091716 

0.058583 

0 . 053583 

0.046394 

0 . 046395 

BETA VALUES FOR 

THE LAYER 

0.000172 

0 . 000035 

-0 . 000063 

0 . 000028 

-0.000022 

0 . 000005 

-0 . 000006 

0 . 000000 

0 . 000000 

- 0.000000 

0.000001 

-0 . 000000 

0.000002 

0.000001 

0.000002 

0 . 000001 

0.000002 

0 . 000000 

0.000001 

0 . 000000 

-0.000003 

-0 . 000000 


ARE : 


3.930280 

3.930289 

3.038528 

3.038530 

1 .761440 

1.761440 

1 ,038440 

1 ,038438 

0.621849 

0.621849 

0.377859 

0.377853 

0.232760 

0.232761 

0.145231 

0 .145231 

0.091716 

0.091716 

0.058583 

0.058583 

0.046394 

0 . 046395 

ARE : 


-0.000020 

-0.000031 

-0 . 000018 

-0.000010 

-0 , 000005 

-0.000007 

0.000001 

-0.000007 

0 . 000003 

-0 . 000002 

0 . 000002 

0 . 000001 

0.000001 

0.000001 

0 .000002 

0 . 000001 

0 ,000001 

0 .000001 

0 . 00000 : 

0.000001 

0 . 00000 ) 

0 . 000001 


3.930271 
3 .038529 
1 . 761452 
1 .038441 
0.621850 
0 . 377860 
0.232762 
0 . 145231 
0.091716 
0 . 058583 
0 .046394 

-0 .000030 
- 0 .000018 
- 0.000010 
- 0.000002 
0 . 000000 
0 ,000001 
0 .000002 
0 ,000002 
0 . 000002 
0.000001 
0 . 000000 



150 


GAMMA VALUES 

for the layer are : 



27,985016 

27.984741 

27.984695 

27.984726 

27.984741 

17.013947 

17.013962 

17.013870 

17.013855 

17.013855 

10.503365 

10 . 503394 

10 . 503348 

10.503335 

10.503335 

6,577075 

6.577109 

6,577087 

6.577077 

6.577081 

4.173593 

4 .173619 

4 . 173608 

4.173602 

4.173599 

2.681605 

2.681616 

2.681611 

2.681610 

2.681606 

1 . 743235 

1 .743243 

1 . 743240 

1 .743237 

1 .743236 

1 .145768 

1 .145776 

1 .145772 

1 .145768 

1 .145767 

0 . 760941 

0.760944 

0.760941 

0.760939 

0 . 760939 

0.510357 

0 .510356 

0 .510354 

0.510353 

0 .510353 

0 . 345494 

0.345491 

0.345490 

0.345491 

0 . 345491 

JACOBIAN VALUES FOR THE LAYER ARE J 



10.487561 

10 . 487476 

10.487505 

10 . 487524 

10.487501 

7. 190097 

7 .190085 

7.190073 

7.190073 

7. 190069 

4 . 301290 

4 . 301294 

4.301282 

4.301279 

4 . 301294 

2 .613405 

2.613412 

2.613409 

2.613404 

2.613409 

1 .611010 

1 .611012 

1 .611010 

1 .61 1009 

1 .611009 

1 . 006612 

1 . 006615 

1 . 006614 

1 . 006613 

1 . 006613 

0,636991 

0.636993 

0.636991 

0,636992 

0.636993 

0.407923 

0 . 407924 

0.407923 

0.407923 

0 . 407923 

0.264179 

0 . 264179 

0.264179 

0.264179 

0.264178 

0.172911 

0. 172911 

0 . 172911 

0. 17291 1 

0. 172910 

0.126605 

0 .126605 

0 .126605 

0.126605 

0 . 126604 
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APPENDIX C 

Computer Output of Spacial Transform 
Coefficients for a Pine Mesh 



2014E-03 0.4000E+00 
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ALPHA VALUES 

FOR THE LAYER 

ARE : 



0.234702 

0 . 234698 

0.234697 

0.234694 

0.234692 

0.222154 

0 . 222152 

0.22214V 

0.222145 

0.222143 

0. 198727 

0. 198724 

0 . 198726 

0.198725 

0.198723 

0. 177777 

0. 177782 

0.177784 

0.177781 

0.177782 

0. 159050 

0. 159053 

0.159054 

0.159053 

0.159055 

0.142305 

0.142304 

0. 142305 

0.142305 

0.142309 

0.127331 

0. 127330 

0 . 127330 

0. 127332 

0.127334 

0.113938 

0 . 1 13938 

0.113933 

0.113941 

0.113942 

0 . 101959 

0.101958 

0.101959 

0 .101959 

0.101961 

0.091242 

0.091242 

0.091242 

0.091243 

0.091244 

0.081655 

0.081655 

0.081655 

0.081656 

0.081657 

0.073079 

0 . 073079 

0.07307V 

0.073080 

0.073079 

0.065406 

0.065407 

0.065405 

0.065407 

0.065406 

0 . 058543 

0.058542 

0.058541 

0.058543 

0.058542 

0.052399 

0.052399 

0.05239? 

0.052400 

0.052399 

0.046903 

0 . 046903 

0.046903 

0.046903 

0.046902 

0.041984 

0.041984 

0.041984 

0.041984 

0.041984 

0.037583 

0.037583 

0.037583 

0.037582 

0.037582 

0.033644 

0 . 033645 

0.033644 

0.033643 

0.033643 

0.030120 

0.030120 

0.030120 

0.030119 

0.030118 
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0.026965 
0.024142 
0.021615 
0.019354 
0.017331 
0.015519 
0.013898 
0 .012446 
0.011147 
0.009983 
0.008941 
0.008009 
0.007173 
0.006425 
0 . 005756 
0.005156 
0 . 004619 
0.004138 
0.003707 
0.003321 
0.003141 


0.026966 
0.024142 
0.021615 
0 .019354 
0.017331 
0.015519 
0.013898 
0.012446 
0.011147 
0.009983 
0.008941 
0.008009 
0.007173 
0.006425 
0.005756 
0.005156 
0.004619 
0.004138 
0.003707 
0.003321 
0.003141 


0.026965 
0.024142 
0.021615 
0.019354 
0.017331 
0 .015520 
0.013898 
0,012447 
0.011147 
0 . 009983 
0.008941 
0 . 008008 
0 .007173 
0 . 006425 
0.005755 
0,005156 
0 .004618 
0.004137 
0.003707 
0.003321 
0.003141 


0.026965 

0.024142 

0.021616 

0.019355 

0.017331 

0.015520 

0.013898 

0.012447 

0.011147 

0.009983 

0.008941 

0.008008 

0.007173 

0.006425 

0.005755 

0.005156 

0.004618 

0.004137 

0.003707 

0.003321 

0.003141 


0.026964 
0.024142 
0 . 021615 
0.019355 
0.017331 
0.015520 
0 .013898 
0.012446 
0.011147 
0.009983 
0.008941 
0 . 008008 
0.007173 
0 . 006425 
0.005755 
0 . 005155 
0.004618 
0.004137 
0.003707 
0.003321 
0.003141 


BETA VALUES FOR THE LAYER ARE : 


0.000075 
- 0.000007 
- 0.000003 
0.000001 
0.000001 
0.000002 
0 . 000003 
0.000003 
0.000002 
0.000002 
0.000002 
0.000002 
0 .000001 
0.000001 
0.000001 
0.000001 
0.000001 
0.000001 
0.000001 
0.000000 


0.000019 
0.000018 
0.000015 
0 . 000010 
0 . 000008 
0 .000007 
0.000005 
0.000004 
0 .000003 
0 . 000002 
0 . 000000 
- 0 . 000000 
- 0.000001 
0 .000000 
-0 . 000000 
- 0 . 000000 
-0 . 000000 
- 0.000001 
- 0.000001 
-0 . 000000 


0.000007 
0.000009 
0.000012 
0.000010 
0 . 000011 
0 . 000011 
0.000009 
0.000007 
0.000005 
0 .000003 
0.000001 
0 .000001 
-0 . OOQGOO 
- 0 . 000001 
- 0 . 000002 
- 0.000002 
- 0 . 000002 
- 0.000002 
- 0 . 000001 
- 0 .000001 


0 . 000002 
0.000005 
0 . 000009 
0.000010 
0.000014 
0.000013 
0.000009 
0.000007 
0 .0 0 0 0 0 £■ 
0.000004 
0.000001 
0 .000000 
- 0 . 0000-00 
- 0.000001 
- 0.000001 
- 0.000002 
- 0.000002 
- 0 .000002 
- 0 .000002 
- 0 . 000000 


0.000006 
0.000001 
0 . 000003 
0 . 000008 
0.000011 
0.000010 
0 . 000007 
0 . 000006 
0 .000005 
0.000003 
0 . 000002 
0.000001 
0.000001 
0.000000 
0 .000000 
- 0.000001 
- 0 . 000000 
- 0.000001 
- 0.000001 
- 0.000001 
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0 . 000000 
0 . 000000 
0.000000 
0 . 000000 
0.000000 
0 . 000000 
0 . 000000 
0 . 000000 
0 . 000000 
0 .000000 
0 . 000000 
0 . 000000 
0 . 000000 
0 . 000000 
-0 . 000000 
- 0 . 000000 
-o. ooooco 
o.ocoooo 
0 . 000000 
0.000000 
- 0.000001 


- 0 . 000000 
0 .000000 
0 . 000000 
0 . 000000 
- 0.000000 
- 0 . 000000 
- 0 . 000000 
- 0.000000 
- 0 . 000001 
- 0 . 000001 
- 0.000001 
- 0 .000000 
- 0 . 000000 
-0 . 000000 
- 0 . 000000 
- 0 . 000000 
- 0.000000 
- 0 . 000000 
0.000000 
- 0.000000 
- 0.000000 


0.000000 
0 . 000001 
0.000001 
0 . 000001 
- 0.000000 
- 0 . 000000 
- 0 . 000000 
- 0 . 000000 
- 0 . 000001 
- 0.000001 
- 0 . 000001 
- 0 .000001 
- 0 . 000001 
- 0 . 000001 
- 0.000000 
- 0 . 000000 
- 0.000000 
0 .000000 
0.000000 
0 . 000000 
0.000000 


0 . 000000 
0 . 000001 
0.000001 
0.000000 
- 0,000000 
- 0 . 000000 
- 0.000000 
- 0.000000 
-0 . 000000 
- 0 . 000000 
- 0.000000 
- 0.000000 
- 0.000001 
- 0.000001 
- 0.000000 
- 0.000000 
- 0.000000 
0 . 000000 
0.000000 
0.000000 
- 0.000000 


- 0.000001 
0 . 000000 
0 . 000001 
0 . 000000 
0 . 000000 
0 , 000000 
- 0 . 000000 
- 0 .000000 
- 0 . 000000 
- 0 .000000 
- 0 .000000 
- 0 . 000000 
- 0 . 000000 
- 0 . 000000 
- 0 .000000 
- 0.000000 
- 0.000000 
0 .000000 
0.000000 
0 . 000000 
- 0 . 000000 



GAMMA VALUES FOR THE LAYER ARE 


1 .982297 

1 . 982210 

1 . 982210 

1.982200 

1 , 982200 

1 . 774610 

1 .774612 

1 .774572 

1 . 774553 

1 . 774546 

1 .588758 

1 . 588765 

1.588745 

1 .588731 

1.588715 

1 .422436 

1 .422443 

1 . 42244 ; 

1 .422431 

1 .422419 

1 . 273581 

1 . 273591 

1,273596 

1.273587 

1 . 273577 

1 . 140352 

1 ,140362 

1 ,140370 

1 .140371 

1 . 140366 

1.021104 

1 .021116 

1 .021130 

1.021133 

1 .021129 

0.914365 

0.914376 

0 . 91438 ° 

0.914397 

0.914396 

0.818816 

0.818828 

0.818841 

0.818851 

0.818850 

0.733284 

0.733292 

0 .733307 

0.733317 

0 .733320 

0.656713 

0 . 656721 

0.656732 

0.656744 

0 . 656748 

0.588163 

0.588168 

0 . 58817 ^ 

0.588190 

0.588195 

0 .526791 

0 . 526794 

0 . 526802 

0.526814 

0 . 526817 

0.471841 

0 .471842 

0.471850 

0 .471860 

0.471864 

0.422640 

0 .422640 

0 . 42264 o 

0.422655 

0.422658 

0.378584 

0 . 378585 

0 . 378589 

0.378596 

0 . 378599 

0.339135 

0.339136 

0.339138 

0.339143 

0.339146 

0.303809 

0.303810 

0 . 303812 

0.303814 

0.303817 

0.272174 

0 . 272175 

0.272176 

0.272177 

0.272179 

0.243843 

0.243844 

0.243844 

0.243845 

0.243846 
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0.218469 
0 .195744 
0, 175389 
0 . 157157 
0. 140826 
0.126197 
0. 113093 
0. 101353 
0.090836 
0.081413 
0 . 072971 
0.065407 
0.058629 
0.052556 
0.047114 
0.042237 
0.037866 
0.033949 
0.030439 
0.027293 
0.024473 


0.218470 
0.195745 
0. 175390 
0.157158 
0 . 140827 
0. 126198 
0.113093 
0 . 101353 
0.090836 
0.081413 
0.072970 
0 . 065406 
0 . 058629 
0.052556 
0.047113 
0.042237 
0 . 037866 
0.033949 
0.030439 
0.027293 
0.024472 


0 .218471 
0 .195745 
0.175391 
0.157159 
0.140828 
0.126199 
0. 113094 
0 . 101354 
0.090836 
0.081413 
0.072970 
0.065406 
0.058628 
0.052555 
0 .0471 13 
0.042236 
0.037865 
0.033948 
0.030438 
0,027292 
0 . 024472 


0.218471 
0.195746 
0.175392 
0.157160 
0.140829 
0.126200 
0.113095 
0. 101355 
0.090837 
0.081414 
0.072971 
0.065406 
0.058628 
0.052554 
0.047112 
0.042235 
0.037865 
0.033948 
0.030438 
0.027291 
0.024472 


0.218472 
0.195747 
0, 175393 
0,157161 
0.140830 
0. 126201 
0 . 113096 
0. 101355 
0.090837 
0.081414 
0.072971 
0.065406 
0 . 058628 
0.052554 
0 . 047112 
0.042235 
0.037864 
0.033948 
0 .030438 
0.027291 
0 .024472 
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JACOBIAN VALUES FOR THE LAYER ARE J 


0.682092 

0.682072 

0.682070 

0.682063 

0.682061 

0.627883 

0.627881 

0.627869 

0.627861 

0.627856 

0.561898 

0.561895 

0.561894 

0.561890 

0.561884 

0.502869 

0.502877 

0.502879 

0.502873 

0.502873 

0.450071 

0 .450076 

0 .450078 

0.450075 

0.450076 

0.402837 

0 .402837 

0 . 40284 I 

0 .402840 

0.402846 

0.360581 

0.360582 

0.360583 

0.360588 

0.360589 

0.322771 

0.322773 

0.322776 

0.322780 

0.322783 

0.288938 

0 .288940 

0 . 288944 

0.288945 

0 . 288948 

0.258662 

0.258664 

0.258667 

0.258670 

0.258671 

0.231568 

0.231570 

0.231572 

0.231575 

0.231577 

0.207322 

0.207323 

0.207325 

0.207327 

0.207327 

0.185621 

0 . 185623 

0 . 185624 

0. 185626 

0 . 185626 

0.166201 

0.166200 

0 .166201 

0.166204 

0 . 166204 

0.148815 

0.148815 

0.148817 

0.148819 

0 .148818 

0.133254 

0 .133254 

0.133254 

0.133256 

0.133256 

0.119324 

0.119324 

0.119325 

0.119325 

0. 119326 

0.106855 

0.106856 

0.106856 

0.106855 

0 . 106856 

0.095693 

0.095694 

0.095693 

0.095692 

0.095692 

0.085700 

0.085700 

0.085700 

0.085699 

0.085698 
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0.076753 

0.068743 

0.061572 

0,055151 

0.049403 

0.044255 

0.039645 

0.035517 

0.031820 

0.028509 

0.025543 

0.022887 

0.020508 

0.018376 

0.016467 

0.014757 

0.013224 

0.011852 

0.010622 

0.009520 

0.008767 


0 . 076754 
0.068743 
0.061572 
0.055151 
0 . 049403 
0.044255 
0 . 039646 
0.035517 
0.031820 
0 . 028509 
0 . 025543 
0.022887 
0.020508 
0.018376 
0.016467 
0.014757 
0.013225 
0.011852 
0.010622 
0.009520 
0.008767 


0.076753 
0.068743 
0.061572 
0.055152 
0 . 049403 
0.044256 
0.039646 
0 . 035518 
0.031821 
0 . 029509 
0 . 025543 
0.022887 
0.020507 
0 .018376 
0.016467 
0 .014756 
0.013224 
0 .011852 
0.010622 
0.009520 
0.008767 


0.076753 
0.068744 
0.061573 
0.055153 
0.049404 
0.044256 
0 . 039646 
0.035518 
0.031820 
0.029509 
0.025543 
0.022887 
0.020507 
0.018376 
0.016466 
0.014756 
0.013224 
0.011851 
0.010622 
0.009520 
0.008767 


0.076752 
0.068743 
0.061573 
0.055153 
0 .049404 
0.044256 
0.039646 
0.035518 
0.031820 
0.028509 
0.025543 
0.022897 
0.020507 
0.018376 
0.016466 
0.014756 
0.013224 
0.011851 
0.010622 
0 . 009520 
0.008767 
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APPENDIX D 


Computer Subroutine for Transformed 
Conduction Equation 
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APPENDIX G 
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Boundary Condition where One Layer’s Boundary 
Undergoes a Change of Phase 
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